Mechanistic PK-PD model of alendronate treatment of postmenopausal osteoporosis predicts bone site-specific response

1 Introduction

Osteoporosis is a disease caused by an imbalance in the remodeling process. Resorption of the old bone, carried out by osteoclasts, predominates over osteoid formation by osteoblasts. This leads to net bone loss and a decrease in the stiffness and strength of bones. Around 9 million osteoporotic fractures occurred worldwide in 2000, with the highest incidence in Europe (Johnell and Kanis, 2006). Over 27.5 million people suffer from this disease in Europe and nearly 3.5 million fractures occur every year, raising healthcare costs to approximately 37 billion euros (Hernlund et al., 2013). Postmenopausal osteoporosis (PMO), resulting from estrogen deficiency, is the most common type of osteoporosis (Eastell et al., 2016).

At present, the most widely used drug for the treatment of PMO is alendronate (Cummings et al., 2020), a bisphosphonate usually administered orally once a week in 70 mg tablets, though other treatment regimens are also prescribed. Alendronate inhibits bone resorption by osteoclasts via two different mechanisms: a reduction in osteoclasts' resorbing capacity and inducing their apoptosis (Halasy-Nagy et al., 2001; Yu et al., 2018; Takagi et al., 2021).

Pharmacokinetics (PK) models intend to explain the absorption, distribution, and elimination of a drug in the body. These PK models constitute the first step to developing pharmacodynamics (PD) models, in which the effect of the drug on the body is analyzed. The combination of both is usually termed as PK-PD models.

Pharmacokinetics of alendronate has been widely studied. Many experimental works have been carried out both on animals and humans. Regarding the latter, there is considerable clinical information about the evolution of the alendronate concentration in plasma (Cocquyt et al., 1999; Chae et al., 2014) and urine excretion (Cocquyt et al., 1999; Kang et al., 2006; Acotto et al., 2012; Chae et al., 2014) during the first hours after alendronate intake. However, to the best of our knowledge, no clinical data have been published regarding the long-term plasma concentration, and only one work has reported on experimental urine excretion in long term (Khan et al., 1997). Despite the available experimental data, only a few mathematical models have been developed to simulate the pharmacokinetics of alendronate or other bisphosphonates (Cremers et al., 2002; Pillai et al., 2004; Sedghizadeh et al., 2013; Chae et al., 2014). These are compartmental models where each compartment represents a part of the body, or rather, a distribution volume where the drug is delivered. These PK models aim at reproducing the temporal evolution of the drug in each compartment. Most of these PK models simply fitted the experimental data from short-term experiments (Pillai et al., 2004; Sedghizadeh et al., 2013; Chae et al., 2014), making their use unreliable for osteoporosis treatments, which normally last for years. Cremers et al. (2002) developed a PK model of pamindronate by fitting the short- and long-term clinical results, but actually they had to use long-term urine excretion data of alendronate (Khan et al., 1997) because the results for pamindronate were not available. To the best of the authors' knowledge, in the case of alendronate treatment, no model could be found in the literature that fits both short- and long-term plasma concentration and/or urine excretion.

Regarding pharmacodynamics, the effects of alendronate on bone density gain (BDG), reduction of fracture risk, and bone turnover markers are well-documented in the literature. However, few mathematical models have been developed to simulate the pharmacodynamics of bisphosphonates and, in particular, alendronate. A large number of these models describe the pharmacodynamics of alendronate in a simplistic way in order to predict the evolution of (systemic) bone turnover markers (Hernández et al., 2002; Pillai et al., 2004). However, they do not consider bone remodeling and its underlying mechanobiological processes. Therefore, these models are unable to predict site-specific bone mineral density (BMD) changes during drug treatment. Recently, Ashrafi et al. (2021) applied a PK-PD model implemented in a bone remodeling model to simulate the behavior of dental implants in patients undergoing bisphosphonate treatment. However, they used ibandronate and implemented the PK model proposed by Pillai et al. (2004), who only used short-term clinical data to adjust their model, which probably makes its use for simulating long-term treatment questionable.

Based on the previous review, the objectives of this work are the following:

• To develop a comprehensive PK model for alendronate valid for short- and long-term treatment.

• To develop a mechanistic PK-PD model of alendronate treatment of PMO to predict site-specific changes in BMD at the hip and lumbar spine.

This article is organized as follows: in Section 2, we provide a detailed description of the different PK and PD models tested in this work and the methods used to fit their constants. The comparison of clinical data with the results provided by the models is reported in Section 3, together with the values of the fitted constants. The results are discussed in detail with respect to the alendronate literature in Section 4.

2 Materials and methods

In order to develop a PK-PD model, it is important to note that the PK and the PD model influence one another, that is, they are not independent. Herein, we first develop the PK model for alendronate with the aim of identifying the most suitable compartment interactions for both short- and long-term responses.¹ At the second stage, we added the PD model which is based on a BCPM of bone remodeling in order to assess changes in BMD at different bone sites. For the case of alendronate treatment of PMO, the addition of the PD will have a greater influence in long term, as the effect of alendronate and its return rate to plasma depends on osteoclasts' action and is more noticeable in long term (months, years) than in short term (hours).

2.1 Pharmacokinetics

The aim of this section was to develop a model to reproduce the pharmacokinetics of alendronate while temporarily neglecting the pharmacodynamics of the drug. Compartmental models will be used to fulfill this objective.

Alendronate is distributed throughout the body via the blood plasma. It directly reaches the plasma in the case of intravenous (IV) administration or through the gastrointestinal tract when delivered orally. It does not appear to be metabolized in humans, and it is cleared from the plasma by deposition onto bone and elimination via renal excretion (Porras et al., 1999). Once in the plasma, alendronate is quickly distributed into the non-calcified tissues of the body, followed by redistribution in bone or renal elimination (Lin et al., 1991; Porras et al., 1999).

Based on that knowledge, a three-compartment model was developed as depicted in Figure 1, where CC stands for the central compartment (plasma), NCT for the non-calcified tissues, and BC for the bone compartment, with the gut compartment added in the case of oral doses. The arrows indicate the flow of alendronate, k_i is the absorption rate constant, and k_el,i is the elimination rate constant, with i = CC, NCT, BC, or urine. The different compartments where alendronate is distributed and the interconnections between them are represented in Figure 1, with the renal excretion represented by an elimination term from the CC (urine excretion). In the following, this model is denoted as a three-compartment model “in parallel” to distinguish it from another model presented later with three compartments and two of them connected “in series.”

FIGURE 1

FIGURE 1. Alendronate three-compartment PK model “in parallel.”

The differential equations governing the temporal evolution of alendronate in the three-compartment PK model “in parallel” are:

$\begin{array}{l}\hfill \begin{array}{cc}\hfill & \frac{dAl{e}_{CC}}{dt}=F\cdot k_{CC}Al{e}_{\mathit{Gut}}+k_{el,BC}Al{e}_{BC}+k_{el,NCT}Al{e}_{\mathit{NCT}}\hfill \\ \hfill & -\left(k_{BC}+k_{\mathit{NCT}}+k_{el,urine}\right)Al{e}_{CC}\hfill \end{array},\end{array}(1)$

$\frac{dAl{e}_{BC}}{dt}=k_{BC}Al{e}_{CC}-k_{el,BC}Al{e}_{BC},(2)$

$\frac{dAl{e}_{\mathit{NCT}}}{dt}=k_{\mathit{NCT}}Al{e}_{CC}-k_{el,NCT}Al{e}_{\mathit{NCT}},(3)$

$\frac{dAl{e}_{\mathit{Gut}}}{dt}=-k_{CC}Al{e}_{\mathit{Gut}},(4)$

$\frac{dAl{e}_{\mathit{Urine}}}{dt}=k_{el,urine}Al{e}_{CC},(5)$

where Ale_CC, Ale_BC, Ale_NCT, Ale_Gut, and Ale_Urine are the amount of alendronate in the central compartment, bone compartment, non-calcified tissues, gut, and urine, respectively. F is the bioavailability. If the dose is administered intravenously, the gut compartment will not appear; therefore, the term F ⋅ k_CCAle_Gut in Eq. 1 must be replaced by the IV dose rate, while Eq. 4 is no longer needed. The concentration in the central compartment was calculated as follows:

$\left[Al{e}_{CC}\right]=\frac{Al{e}_{CC}}{V_c},(6)$

where V_c is the volume of distribution of the central compartment.

Other authors have proposed models with different numbers of compartments and different compartmental interactions. In order to find the most suitable PK model which minimizes the error between experimental data and model results, we investigated four additional models described below:

(1) One-compartment model with two elimination mechanisms, proposed by Chae et al. (2014) (see Figure 2A). This model has an important drawback. It is known from the literature that alendronate is only excreted by urine (Porras et al., 1999), but the model consists of a CC and two exits: one to urine and another to a non-specified organ (in the original model, the outflow rates were termed k_urine and k_non−urine). The latter was renamed here as k_BC in order to have an analogous terminology in all models. To overcome this drawback, the following two-compartment model was proposed.

(2) Two-compartment model. The only difference with the one-compartment model is that one of the exits from the CC flows into the BC, and there is a return flow from BC to CC (see Figure 2B).

(3) Three-compartment model “in series.” This model adds to the previous model, a compartment connected to the BC, as seen in Figure 2C. The model was called “in series” to distinguish it from the three-compartment model “in parallel.” This model considers the different availability of the alendronate deposited in bone. To this end, it distinguishes the drug that has been deposited near the bone matrix–marrow interface from the alendronate that was buried deeper into the bone matrix. The former, termed active alendronate (and contained in the BC), is more accessible for osteoclasts to resorb it since resorption occurs mainly on the bone matrix–marrow interface. Thus, it can affect the osteoclastic activity through endocytosis. On the contrary, the latter is assumed to be inaccessible to resorption and was termed inactive alendronate (contained in the inactive compartment, IC) (Porras et al., 1999). Cremers et al. (2002) used a similar model for pamindronate.

(4) Four-compartment model. It was proposed by Porras et al. (1999), although they did not test it. It can be considered the combination of the three-compartment models “in parallel” and “in series,” where the biological backgrounds of both apply (see Figure 2D).

FIGURE 2

FIGURE 2. Different types of PK models for alendronate considering various numbers of compartments and different compartmental interactions. (A) One-compartment model with two elimination mechanisms. (B) Two-compartment model. (C) Three-compartment model “in series.” (D) Four-compartment model.

The differential equations governing the previous models are not included in the main text but can be consulted in the Supplementary Material. In accordance with the literature (Lin et al., 1992; Cocquyt et al., 1999), no saturation has been considered for the bone compartment in any of the PK models, though the concept of saturation will be introduced later when the fully mechanistic PK-PD model is developed.

The information found in the literature about the pharmacokinetics of alendronate consists of clinical data of the serum concentration in short term and urine excretion amount in short and long terms. Little is known about urinary excretion in long term, and only the experimental data obtained by Khan et al. (1997) could be found to validate these specific results. The experimental data from Chae et al. (2014) were used to adjust the serum concentration in short term, while those obtained by Kang et al. (2006) were used for the short-term urinary excretion. Other authors measured similar values for both types of results (Cocquyt et al., 1999; Acotto et al., 2012). It is important to note that the results from Chae et al. and Kang et al. correspond to a single oral administration (70 mg of alendronate) to healthy volunteers, while the results from Khan et al. correspond to IV administration (7.5 mg of alendronate infused over 12 h every day for 4 days, being a total dose of 30 mg) to postmenopausal women. These dosing regimens were simulated with the PK models, and their constants were fitted to reproduce the clinical results (short-term serum concentration and short-term and long-term urinary excretion). The goodness of fit of the three curves was jointly evaluated to identify the most suitable PK model. A genetic algorithm was used to minimize the difference between the curves generated by the model and the clinical results. Since the three sets of clinical results were given in different units, the root mean squared error of each curve was normalized with the mean of each experimental curve in order to express the errors as a coefficient of variation so that they could be summed up. Finally, this sum was divided by 3 and expressed in percentage to yield E_PK. This error measure represents the average error over all three experiments and was the objective function minimized with the genetic algorithm:

$\begin{array}{cc}\hfill E_{PK}=& \left(\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Chae}}}{\left(S_{\mathit{model}}\left(t_i\right)-S_{\mathit{Chae}}\left(t_i\right)\right)}^2}{N_{\mathit{Chae}}}}}{{\mu }_{\mathit{Chae}}}+\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Kang}}}{\left(s{U}_{\mathit{model}}\left(t_i\right)-s{U}_{\mathit{Kang}}\left(t_i\right)\right)}^2}{N_{\mathit{Kang}}}}}{{\mu }_{\mathit{Kang}}}+\right.\hfill \\ \hfill & +\left.\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Khan}}}{\left(l{U}_{\mathit{model}}\left(t_i\right)-l{U}_{\mathit{Khan}}\left(t_i\right)\right)}^2}{N_{\mathit{Khan}}}}}{{\mu }_{\mathit{Khan}}}\right)\cdot \frac{100}{3}\hfill \end{array}(7)$

where N_j is the number of time points and μ_j is the mean of experimental values over the respective time points for each set of clinical data (j = Chae, Kang, Khan). The variables S_j(t_i), sU_j(t_i), and lU_j(t_i) represent the serum concentration, the short-term urinary excretion, and long-term urinary excretion at the time point t_i of the clinical result j, respectively. Finally, the subscript j = model stands for the prediction of the PK model.

2.2 Pharmacodynamics

Once a suitable PK model is identified, the effect of alendronate on bone turnover was modeled. To address this aim, a previously published mathematical BCPM describing the bone cell interactions was used (Martin et al., 2019). This model considers catabolic (RANK–RANKL–OPG) and anabolic (Wnt–Scl– LRP5/6) signaling pathways, together with the action of parathyroid hormone (PTH), nitric oxide (NO), transforming growth factor beta (TGF-β), and mechanobiological feedback on bone cells. The effect of bone mineralization was added following Martínez-Reina and Pivonka (2019) and Martínez-Reina et al. (2008). The accumulation and repair of microstructural damage was also taken into account as in Martínez-Reina et al. (2021).

The bone cell types, whose concentrations are the state variables of the model, are osteoblast precursor cells (Ob_p), active osteoblasts (Ob_a), osteoclast precursor cells (Oc_p), active osteoclasts (Oc_a), and osteocytes (Ot). The cell pools of uncommitted osteoblasts (Ob_u) and osteoclasts (Oc_u) are assumed constant as in Martin et al. (2019).

$\begin{array}{l}\hfill \begin{array}{cc}\hfill & \frac{d{\mathrm{O}\mathrm{b}}_{\text{p}}}{dt}=D_{{\mathrm{O}\mathrm{b}}_{\text{u}}}\cdot {\mathrm{O}\mathrm{b}}_{\text{u}}\cdot {\mathrm{\pi }}_{\text{act,Obu}}^{\text{TGF}-\mathrm{\beta }}-D_{{\mathrm{O}\mathrm{b}}_{\text{p}}}\cdot {\mathrm{O}\mathrm{b}}_{\text{p}}\cdot {\mathrm{\pi }}_{\text{rep},{\mathrm{O}\mathrm{b}}_{\text{p}}}^{\text{TGF}-\mathrm{\beta }}\hfill \\ \hfill & +P_{{\mathrm{O}\mathrm{b}}_{\text{p}}}\cdot {\mathrm{O}\mathrm{b}}_{\text{p}}\cdot {\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{b}}_{\text{p}}}^{\text{Wnt}}\hfill \end{array},\end{array}(8)$

$\frac{d{\mathrm{O}\mathrm{b}}_{\text{a}}}{dt}=D_{{\mathrm{O}\mathrm{b}}_{\text{p}}}\cdot {\mathrm{O}\mathrm{b}}_{\text{p}}\cdot {\mathrm{\pi }}_{\text{rep},{\mathrm{O}\mathrm{b}}_{\text{p}}}^{\text{TGF}-\mathrm{\beta }}-{\mathrm{\Delta }}_{{\mathrm{O}\mathrm{b}}_{\text{a}}}\cdot {\mathrm{O}\mathrm{b}}_{\text{a}},(9)$

$\frac{d{\mathrm{O}\mathrm{c}}_{\text{p}}}{dt}=D_{{\mathrm{O}\mathrm{c}}_{\text{u}}}\cdot {\mathrm{O}\mathrm{c}}_{\text{u}}\cdot {\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{u}}}^{\text{RANKL}}-D_{{\mathrm{O}\mathrm{c}}_{\text{p}}}\cdot {\mathrm{O}\mathrm{c}}_{\text{p}}\cdot {\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{p}}}^{\text{RANKL}},(10)$

$\frac{d{\mathrm{O}\mathrm{c}}_{\text{a}}}{dt}=D_{{\mathrm{O}\mathrm{c}}_{\text{p}}}\cdot {\mathrm{O}\mathrm{c}}_{\text{p}}\cdot {\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{p}}}^{\text{RANKL}}-A_{{\mathrm{O}\mathrm{c}}_{\text{a}}}\cdot {\mathrm{O}\mathrm{c}}_{\text{a}}\cdot {\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{p}}}^{\text{TGF}-\mathrm{\beta }},(11)$

$\frac{d\mathrm{O}\mathrm{t}}{dt}=\eta \frac{d{f}_{bm}}{dt},(12)$

where $D_{{\mathrm{O}\mathrm{b}}_{\text{u}}}$, $D_{{\mathrm{O}\mathrm{b}}_{\text{p}}}$, $D_{{\mathrm{O}\mathrm{c}}_{\text{u}}}$, and $D_{{\mathrm{O}\mathrm{c}}_{\text{p}}}$ are the differentiation rates of Ob_u, Ob_p, Oc_u, and Oc_p, respectively; $A_{{\mathrm{O}\mathrm{c}}_{\text{a}}}$ is the apoptosis rate of Oc_a, and ${\mathrm{\Delta }}_{{\mathrm{O}\mathrm{b}}_{\text{a}}}$ is the rate of clearance of active osteoblasts through apoptosis or differentiation into osteocytes. The variables ${\mathrm{\pi }}_{\text{act,Obu}}^{\text{TGF}-\mathrm{\beta }}$, ${\mathrm{\pi }}_{\text{rep},{\mathrm{O}\mathrm{b}}_{\text{p}}}^{\text{TGF}-\mathrm{\beta }}$, and ${\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{p}}}^{\text{TGF}-\mathrm{\beta }}$ represent activator and repressor functions related to the binding of TGF-β to its receptor. Similarly, ${\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{u}}}^{\text{RANKL}}$ and ${\mathrm{\pi }}_{\text{act},{\mathrm{O}\mathrm{c}}_{\text{p}}}^{\text{RANKL}}$ are the activator functions related to the RANK–RANKL binding. Finally, $P_{{\mathrm{O}\mathrm{b}}_{\text{p}}}$ is the proliferation rate of Ob_p, a process which is mediated by the Wnt signaling pathway through the activator function ${\mathrm{\pi }}_{act,O{b}_p}^{\mathit{Wnt}}$. These functions, as well as the remaining equations needed to complete the model, are not presented here for brevity but are completely described in the Supplementary Material. A schematic figure of the mechanistic PK-PD model is presented in Figure 3. The values of the model constants are given in Supplementary Table S1.

FIGURE 3

FIGURE 3. Schematic representation of the mechanistic PK-PD model: bone cell differentiation stages along with biochemical and biomechanical interactions are presented. The mineralization of osteoid is shown in orange. Alendronate doses lead to the distribution of the drug into the bone compartment where it interacts with the osteoclasts.

Finally, Eq. 12 establishes that the population of osteocytes varies as the bone matrix fraction f_bm, if the density of osteocytes trapped within the bone matrix, η, is assumed constant as carried out in Martin et al. (2019). Bone matrix fraction is defined as the volume of the bone matrix, V_b, per total volume of the bone sample (i.e., the representative volume element, V_RVE), expressed as a percentage:

$f_{bm}\left(\%\right)=\frac{V_b}{V_{\mathit{RV\; E}}}\cdot 100.(13)$

Its evolution is obtained through the balance between resorbed and formed tissue:

$\frac{d{f}_{bm}}{dt}=-k_{\mathit{res}}\cdot {\mathrm{O}\mathrm{c}}_{\text{a}}+k_{\mathit{form}}\cdot {\mathrm{O}\mathrm{b}}_{\text{a}},(14)$

where k_res and k_form are the rates of bone resorption and osteoid formation, respectively (see Supplementary Table S1).

This BCPM was integrated into the PK model by including it in the bone compartment (BC, see Figure 1). All the absorption rate constants k_i and the elimination rate constants k_el,i fitted for the PK model (see Figure 1) were assumed constant, except for k_el,BC, which is variable and depends on the osteoclast-mediated release of alendronate from the bone matrix. Note that the latter variable depends on bone turnover (Fisher et al., 2000; Yu et al., 2018), whereas in the PK model of Section 2.1, k_el,BC was constant. Now, in the PK-PD model, k_el,BC is given by:

$k_{el,BC}=k_{\mathit{res}}O{c}_a.(15)$

We have assumed that it is the concentration of alendronate in the bone compartment, BC, that controls the action of the drug on the bone remodeling process:

$\left[Al{e}_{BC}\right]=\frac{Al{e}_{BC}}{V_{\mathit{Bone}}},(16)$

where V_Bone is the total volume of bone tissue in the body, 2.23 ⋅ 10⁻³m³, adapted from Valentín (2002) for a 60 kg adult female. The amount of alendronate accumulated within the bone increases with each weekly dose; therefore, if the effect of the drug was proportional to that amount, such an effect would be increasingly pronounced, but this is not observed in the clinical trials (Rizzoli et al., 2002), where the increase of BMD is high during the first months after the beginning of the treatment and slows down in the subsequent months, both in the lumbar spine and the hip. This behavior could be explained by several facts highlighted in the literature:

• Alendronate is preferentially deposited in areas undergoing active resorption (Porras et al., 1999).

• 80% of normal bone turnover occurs in trabecular bone and 20% in cortical bone (Fleisch, 2000).

• A considerably higher amount of alendronate is deposited in trabecular bone than in cortical bone (Lin et al., 1991; Khan et al., 1997).

• The volume of cortical bone is much higher than that of trabecular bone (Fleisch, 2000; Valentín, 2002).

Therefore, a higher proportion of alendronate would be initially deposited in trabecular bone, which explains the fast increase of BMD in the lumbar spine and the hip observed in the clinical results. However, due to the fact that the volume of trabecular bone is much lower than that of cortical bone, in a long treatment, it is expected that the proportion of alendronate that reaches the former will be reduced over time in favor of the latter due to the saturation of the tissue. This would explain the long-term stabilization of bone mass gain observed clinically.

In view of the clinical results, we hypothesized that the BC can be divided into two parts: one part termed active, which is the closest to the bone marrow interface and where the retrieval of alendronate is immediate through bone resorption; the second part termed inactive corresponds to the innermost regions of bone, where alendronate is buried and, to some extent, inaccessible to bone resorption (Sato et al., 1991; Porras et al., 1999). An updated compartmental model is shown in Figure 4. The active subcompartment would predominate in trabecular bone, where most of the tissue is superficial, and its proportion would decrease with porosity as the tissue becomes more compact.

FIGURE 4

FIGURE 4. Three-compartment model “in parallel”, adapted with the division of the BC compartment into an active (A) and an inactive (I) subcompartment.

We need a variable to measure this compacity, that is, the proportion of superficial tissue. The specific surface, S_v, is defined as the area of bone matrix–marrow interface, S_i, per total volume of the bone sample, V_T. The following correlation between S_v and vascular porosity, f_vas, was given by Martin (1984) and used by Pivonka et al. (2013):

$S_v=\frac{S_i}{V_T}=32.3\cdot f_{\mathit{vas}}-93.9\cdot f_{\mathit{vas}}^2+134\cdot f_{\mathit{vas}}^3-101\cdot f_{\mathit{vas}}^4+28.8\cdot f_{\mathit{vas}}^5,(17)$

where the vascular porosity is $f_{\mathit{vas}}=1-\frac{f_{bm}}{100}$ and S_v is expressed in mm²/mm³, being this variable equal to 0 both for f_vas = 0 and f_vas = 1 when no free surface exists. In order to measure the compacity of the tissue, we would need to express that specific surface per bone matrix volume instead of per total volume:

$\frac{S_i}{V_b}=\frac{\frac{S_i}{V_T}}{\frac{V_b}{V_T}}=\frac{S_v}{f_{bm}}.(18)$

This quotient is not defined for f_bm = 0, but Eq. 18 is not needed for very low values of f_bm as we can assume that all the tissue is superficial and accessible to bone resorption in a very porous bone. More precisely, we have established f_bm0 = 5% as the minimum bone matrix below which all the tissue is considered superficial. We have normalized the quotient in Eq. 18 and used the following expression for the active alendronate, that is, the amount of alendronate contained in the active subcompartment:

$Al{e}_{BC,act}=Al{e}_{BC}\cdot \frac{S_v}{f_{bm}}\cdot \frac{f_{\mathit{bm0}}}{S_{v0}},(19)$

where S_v0 is the specific area corresponding to the reference value f_bm0. The concentration of active alendronate can be assessed as in Eq. 16, that is, $[Al{e}_{BC,act}]=\frac{Al{e}_{BC,act}}{V_{\mathit{bone}}}$. The function $\frac{S_v}{f_{bm}}\cdot \frac{f_{\mathit{bm0}}}{S_{v0}}$ that controls the division between the active and inactive parts of the BC is plotted against f_bm in Figure 5. If the tissue has no pores (f_bm = 100%), all the BC is inactive as there is no surface where osteoclasts can resorb bone; thus, there is no release of alendronate from the bone matrix. As porosity increases (f_bm decreases), the ratio between the free surface and bone volume rises, as does the exposure of the drug (the active subcompartment becomes predominant). In contrast, if more tissue is formed than resorbed, f_bm increases and the inactive part becomes predominant as the alendronate that was on the surface is buried into the bone matrix and is thus less accessible to osteoclasts.

FIGURE 5

FIGURE 5. Function that regulates the division between the active (A) and inactive (I) parts of BC.

So, Eqs 1 and 2 are rewritten as follows by changing the second terms on the right-hand side of both:

$\begin{array}{l}\hfill \begin{array}{cc}\hfill & \frac{dAl{e}_{CC}}{dt}=F\cdot k_{CC}Al{e}_{\mathit{Gut}}+k_{el,BC}\left[Al{e}_{BC,act}\right]\frac{V_{\mathit{Bone}}}{\frac{f_{bm}^{\mathit{aver}}}{100}}+k_{el,NCT}Al{e}_{\mathit{NCT}}\hfill \\ \hfill & -\left(k_{BC}+k_{\mathit{NCT}}+k_{el,urine}\right)Al{e}_{CC}\hfill \end{array},\end{array}(20)$

$\frac{dAl{e}_{BC}}{dt}=k_{BC}Al{e}_{CC}-k_{el,BC}\left[Al{e}_{BC,act}\right]\frac{V_{\mathit{Bone}}}{\frac{f_{bm}^{\mathit{aver}}}{100}}.(21)$

Since [Ale_BC,act] is the local concentration, k_el,BC[Ale_BC,act] represents the local amount of alendronate per unit volume or the RVE, V_RVE, released from the bone matrix through resorption. As the flux from BC to CC is systemic, one must consider the contribution of all the V_RVE in the skeleton, which is done in a simplified way through the factor $V_{\mathit{Bone}}/(f_{bm}^{\mathit{aver}}/100)$. Therefore, $f_{bm}^{\mathit{aver}}/100$ represents an average bone volume fraction of the skeleton that allows expressing the amount of alendronate per unit volume of bone tissue (recall Eq. 13), which is then multiplied by the total volume of bone tissue in the skeleton, V_Bone. Eqs 20 and 21 imply that the systemic flux from BC to CC is proportional to the flux from the local RVE to the CC; if bone turnover changes locally for any reason, it will change accordingly in the whole skeleton. The adopted value $f_{bm}^{\mathit{aver}}=43.7\%$ was estimated by assuming an average f_bm = 93% for cortical bone (Cardoso et al., 2013) and f_bm = 14% for trabecular one (Ulrich et al., 1999). Therefore, 80% of the bone tissue volume is cortical, and the rest is trabecular (Valentín, 2002).

The active subcompartment is predominant in trabecular bone, whose volume is considerably lower than that of cortical bone. Thus, it is plausible to assume that this subcompartment will become saturated relatively soon, whereas the total alendronate in the BC will not, as indicated in the previous subsection because there is a large volume of cortical bone in the body able to admit high doses of the drug. It could also become saturated in cases of very long treatments, and then the excess of alendronate would be eliminated via urine, but the authors found no information on urine excretion in long treatments. The active alendronate has also been saturated as a function of S_v/f_bm, being its maximum concentration:

${\left[Al{e}_{BC,act}\right]}^{max}=f\cdot \frac{S_v}{f_{bm}}\cdot \frac{f_{\mathit{bm0}}}{S_{v0}},(22)$

with f as a constant to be fitted. Thus, as alendronate enters the BC (through Eq. 21), it is divided between the active and inactive subcompartments following Eq. 19 ² until $[Al{e}_{BC,act}]={[Al{e}_{BC,act}]}^{\max }$ when the active subcompartment becomes saturated and no more drug is allowed in it.

Alendronate has the following two effects on bone cells. On the one hand, it causes the disappearance of clear zones and ruffled borders, disrupting the cytoskeleton of osteoclasts by inhibiting farnesyl pyrophosphate (FPP) synthase, which leads to these structural changes and loss of function (Halasy-Nagy et al., 2001; Takagi et al., 2021). In other words, it limits the resorbing capacity of osteoclasts, which is measured in the model by k_res (see Eq. 14). Consequently, this parameter is reduced as follows:

$k_{\mathit{res}}=k_{\mathit{res,nom}}\left(1-{\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}\cdot k_{el,BC}\cdot \left[Al{e}_{BC,act}\right]\right),(23)$

where k_res,nom is the nominal value of k_res, ${\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}$ is a constant quantifying the effect of alendronate on the resorbing capacity of osteoclasts, and k_el,BC ⋅ [Ale_BC,act] is the amount of alendronate per unit volume released from the bone matrix through resorption, that is, the concentration of alendronate that affects the surrounding osteoclasts. If Eq. 15 is replaced in Eq. 23:

$k_{\mathit{res}}=k_{\mathit{res,nom}}\left(1-{\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}\cdot k_{\mathit{res}}\cdot O{c}_a\cdot \left[Al{e}_{BC,act}\right]\right),(24)$

from which k_res can be worked out:

$k_{\mathit{res}}=\frac{k_{\mathit{res,nom}}}{1+{\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}\cdot k_{\mathit{res,nom}}\cdot O{c}_a\cdot \left[Al{e}_{BC,act}\right]}.(25)$

The second effect caused by alendronate is the inhibition of the FPP synthase in the mevalonate pathway and the reduction of protein prenylation, an essential post-translational lipid modification required for the function of numerous proteins, thereby inducing apoptosis in osteoclasts (Sato et al., 1991; Halasy-Nagy et al., 2001; Yu et al., 2018; Takagi et al., 2021). To account for this effect, the apoptosis rate was increased in the model as follows:

$A_{O{c}_a}=A_{O{c}_a,nom}\left(1+{\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}\cdot k_{el,BC}\cdot \left[Al{e}_{BC,act}\right]\right),(26)$

where ${\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}$ quantifies the effect of alendronate on the apoptosis of osteoclasts. If Eq. 15 is replaced in Eq. 26:

$A_{O{c}_a}=A_{O{c}_a,nom}\left(1+{\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}\cdot k_{\mathit{res}}O{c}_a\cdot \left[Al{e}_{BC,act}\right]\right).(27)$

Therefore, the constants to be fitted in the PK-PD model are k_CC, k_NCT, k_BC, k_el,NCT, k_el,urine, V_c, and F, the purely PK constants of the three-compartment model “in parallel”, which was the model eventually selected, and ${\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}$, ${\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}$, and f, the PD constants that affect the PK model through k_el,BC and Eq. 15. It should be noted that due to this coupling through k_el,BC, the PK model constants needed to be readjusted. The same clinical data used for the fitting of the PK model were used to validate the PK-PD model: Khan et al. (1997) for the urinary excretion in the long term, Chae et al. (2014) for the serum concentration in the short term, and Kang et al. (2006) for the short-term urinary excretion. Furthermore, the effect of alendronate on the bone response, included in the PD model, had to be validated with bone turnover biomarkers. In this case, only local markers could be used, as the present model was applied locally, at a specific bone location that is only characterized by its BMD. For this reason, the bone density gain (BDG), measured with respect to baseline (the beginning of the treatment when bone apparent density was ρ₀), was defined as:

$\mathrm{B}\mathrm{D}\mathrm{G}\left(\%\right)=\frac{\rho \left(t\right)-{\rho }_0}{{\rho }_0}\cdot 100.(28)$

This BDG was compared to clinical data taken from Rizzoli et al. (2002), who measured the BDG in the lumbar spine (for which f_bm = 15% was assumed) and the hip (f_bm = 25%) in postmenopausal women with osteoporosis treated with a once-weekly dose of 70 mg of alendronate during 2 years. PMO was simulated in the in silico model as in Martínez-Reina et al. (2022) by introducing a disease-related increase in RANKL production over time, $P_{\mathit{RANKL}}^{\mathit{PMO}}$, through the following sigmoidal function:

$P_{\mathit{RANKL}}^{\mathit{PMO}}=P_{\mathit{RANKL}}^{\mathit{PMO,max}}\frac{{\left(t-t_{\mathit{onset}}\right)}^{\gamma }}{{\left(t-t_{\mathit{onset}}\right)}^{\gamma }+{\delta }_{\mathit{PMO}}^{\gamma }},(29)$

where $P_{\mathit{RANKL}}^{\mathit{PMO,max}}$ is the maximum (long term) RANKL production rate due to PMO, t_onset is the onset of RANKL increase, τ_PMO is a time constant that establishes when the 50% of $P_{\mathit{RANKL}}^{\mathit{PMO,max}}$ is reached, and δ_PMO is a parameter which controls how steep the increase of RANKL is (i.e., the resemblance to a step function). The time between the onset of RANKL increase and the beginning of the treatment was set to 19 years, which is the average time elapsed since menopause in the patients studied by Rizzoli et al. (2002).

The PK-PD model was fitted using a genetic algorithm to minimize the difference between the curves generated by the model and the clinical results. The definition of the error is similar to the one used in the fitting of the PK model (see Eq. 7) but with two differences: 1) the BDG results were also compared, in this case, to those obtained by Rizzoli et al. (2002); 2) it was necessary to distinguish between two types of bone and evaluate them separately (f_bm = 15%, assumed for the lumbar spine, and f_bm = 25% assumed for the hip) as Rizzoli et al. (2002) did in their clinical measurements of BDG. This distinction had to be carried out also for the serum concentration and short- and long-term urinary excretion since the corresponding in silico results depend on f_bm (see Eqs 15, 25, and 27). In fact, Eq. 15 couples the PK and PD models and makes the prediction of the serum and urine profiles also dependent on f_bm. Thus, the error defined as the objective function and minimized with the genetic algorithm was

$\begin{array}{cc}\hfill & E=\sum _{k=1}^2\left(\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Chae}}}{\left(S_{\mathit{model,k}}\left(t_i\right)-S_{\mathit{Chae}}\left(t_i\right)\right)}^2}{N_{\mathit{Chae}}}}}{{\mu }_{\mathit{Chae}}}+\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Kang}}}{\left(s{U}_{\mathit{model,k}}\left(t_i\right)-s{U}_{\mathit{Kang}}\left(t_i\right)\right)}^2}{N_{\mathit{Kang}}}}}{{\mu }_{\mathit{Kang}}}+\right.\hfill \\ \hfill & +\left.\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Khan}}}{\left(l{U}_{\mathit{model,k}}\left(t_i\right)-l{U}_{\mathit{Khan}}\left(t_i\right)\right)}^2}{N_{\mathit{Khan}}}}}{{\mu }_{\mathit{Khan}}}+\frac{\sqrt{\frac{\sum _{i=1}^{N_{\mathit{Riz,k}}}{\left(BD{G}_{\mathit{model,k}}\left(t_i\right)-BD{G}_{\mathit{Riz,k}}\left(t_i\right)\right)}^2}{N_{\mathit{Riz,k}}}}}{{\mu }_{\mathit{Riz,k}}}\right)\cdot \frac{100}{8}\hfill \end{array}(30)$

where k = 1, 2 correspond to f_bm = 15% and f_bm = 25%, respectively, N_j is the number of time points, and μ_j is the mean of those points for each set of clinical data (j = Chae, Kang, Khan, Riz). The variables S_j(t_i), sU_j(t_i), lU_j(t_i), and BDG_j(t_i) represent the serum concentration, the short- and long-term urinary excretion, and the BDG at the time point t_i, respectively. Finally, the subscript j = model stands for the prediction of the PK-PD model for each type of bone k.

A sensitivity analysis was conducted to evaluate the uncertainty of the model predictions. The fitted constants of the PK-PD model were varied by ± 10%, one at a time, keeping the rest constant.

In order to demonstrate the implications of different assumptions made for the development of the PK-PD model and more specifically, in the distribution of alendronate and its effect on bone turnover, six scenarios have been tested with six different models:

(1) PK-PD1. This is the proposed model for which alendronate is assumed to downregulate the resorption capacity of osteoclasts (through Eq. 25) and upregulate its apoptosis rate (through Eq. 27) as confirmed by different authors (Halasy-Nagy et al., 2001; Takagi et al., 2021). It is also hypothesized that the alendronate in BC is divided into active and inactive subcompartments and that the active one can be saturated.

(2) PK-PD2. In this model, we assume that alendronate only affects the resorption capacity but not the apoptosis rate ($A_{O{c}_a}=A_{O{c}_a,nom}$ in Eq. 27).

(3) PK-PD3. Herein, we assume that alendronate only affects the apoptosis rate but not the resorption capacity (k_res = k_res,nom in Eq. 25).

(4) PK-PD4. In this model, we assume that the active subcompartment does not exhibit saturation. Therefore, Eq. 22 is removed from this model.

(5) PK-PD5. Herein, we assume that all the alendronate stored in the BC (and not just the active one) is accessible to osteoclasts and therefore influences the PD model. This is equivalent to disregarding the distinction between active and inactive subcompartments and replacing [Ale_BC,act] by [Ale_BC] in Eqs 25 and 27.

(6) PK-PD6. As in the previous model PK-PD5, we assume here that all the alendronate ([Ale_BC]) influences the PD model but including saturation similar to models PK-PD1, PK-PD2, and PK-PD3 for the active subcompartment. As no distinction is made here between active and inactive, when the total alendronate reaches the maximum concentration, no further alendronate can accumulate in the BC. This is in contrast to models PK-PD1, PK-PD2, and PK-PD3, where alendronate can still enter the inactive subcompartment when the active one is saturated.

The constants of the six models were fitted separately by minimizing the error E defined in Eq. 30 with a genetic algorithm. The fitting of models PK-PD2 and PK-PD3 was made by enforcing ${\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}=0$ and ${\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}=0$, respectively.

Finally, once the constants were fitted to show that the PK-PD1 model was the most accurate, this was used to predict the BDG in postmenopausal women with osteoporosis treated with a daily dose of 10 mg of alendronate, and the results were compared with the clinical data from Rizzoli et al. (2002).

3 Results

The results of the fitting used to determine the most suitable PK model are presented first. Subsequently, the results of integrating this PK model with the various PD models are shown.

3.1 Pharmacokinetics

In Table 1, the goodness of fit of each PK model (through error E_PK, see Eq. 7) is presented, together with the values of the fitted constants. In Figure 6, the fitting curves of each PK model to the clinical results of plasma concentration and short- and long-term urine excretion are shown.

TABLE 1

TABLE 1. Summary of the constants and errors of the alendronate PK models. 1C: one-compartment model with two exits; 2C: two-compartment model; 3C series: three-compartment model “in series”; 3C parallel: three-compartment model “in parallel”; 4C: four-compartment model.

FIGURE 6

FIGURE 6. PK model results of alendronate: (A) plasma concentration versus time; (B) short-term urinary excretion vs. time; and (C) long-term urinary excretion vs. time. 1C: one-compartment model with two exits; 2C: two-compartment model; 3C series: three-compartment model “in series”; 3C parallel: three-compartment model “in parallel”; 4C: four-compartment model.

As can be seen in Figure 6, all the PK models are able to reproduce with good accuracy the temporal evolution of the concentration of alendronate in the plasma (Figure 6A) and the short-term urinary excretion (Figure 6B), with negligible differences between them. However, the one-compartment model with two exits and the two-compartment model are unable to reproduce urinary excretion in the mid and long terms (see Figure 6C). In the case of the one-compartment model with two exits (blue solid line), the alendronate is eliminated from the central compartment via urine, and once it is exhausted, no more can be eliminated, thus leading to a step-like curve. The two-compartment model (black dashed line) shows a significant improvement with respect to the one-compartment model: it yields a curve with a non-zero slope in the mid and long terms, indicating that alendronate continues to be eliminated in those periods. However, that slope does not fit the real one. An option for the fitting procedure was tried in this case: enforcing the long-term slope, but then the model completely overestimated the amount of alendronate excreted in the midterm (30–100 days). The rest of the models (three-compartment “in series” and “in parallel” and four-compartment) produce almost identical results, which are all in excellent agreement with the experimental data.

Table 1 shows that the one-compartment model with two exits and the two-compartment model produced errors considerably higher than the rest of the models, which in turn have similar errors. Therefore, it is important to highlight the two qualitative leaps made when moving from 1 to 2 and from 2 to 3 compartments. In view of these results, it seems that at least three compartments are needed to simulate correctly the alendronate elimination rate over time.

The three-compartment model “in parallel” was selected as the most suitable PK model, because it is simple and able to fit the three experimental curves with good accuracy. The selected model was only slightly outperformed by the four-compartment model, but the difference was virtually negligible, and the latter model is more complex as it depends on two more constants. Moreover, the layout of the three-compartment model “in parallel” is biologically based and allows a direct interpretation of all the compartments and their mutual interaction.

3.2 Pharmacodynamics

In Table 2, the final constants fitted for the PK-PD model are presented. If the constants of the PK model (without PD) are compared with those provided in Table 1 for the three-compartment PK model “in parallel,” it can be seen that they are very similar, with slight variations due to coupling with the PD model.

TABLE 2

TABLE 2. Summary of parameter values for alendronate PK models. Fitted for the complete model (PK-PD1), for the model that does not take into account the influence of alendronate on the apoptosis of osteoclasts (PK-PD2), for the model that does not consider its influence on the resorbing capacity of osteoclasts (PK-PD3), for the model in which the active subcompartment does not exhibit saturation (PK-PD4), for the model that ignores the distinction between the active and inactive subcompartments (PK-PD5), and for the model that ignores that distinction but implements saturation in the BC (PK-PD6). E is the error as defined in Eq. 30.

In Figure 7, the plasma concentration and the short- and long-term urine excretion curves obtained with the model PK-PD1 are shown along with the experimental data for f_bm = 15% (results for f_bm = 25% are omitted because the curves virtually overlap with those of f_bm = 15%, and only in the long-term urine excretion, a slight difference can be noted in the final slope). It can be seen that the model PK-PD1 is able to reproduce, reasonably well, all the experimental curves.

FIGURE 7

FIGURE 7. PK-PD1 model results of alendronate: (A) plasma concentration versus time; (B) short-term urinary excretion vs. time; and (C) long-term urinary excretion vs. time.

In Figure 8, the BDG predicted by the model PK-PD1 is plotted against the experimental results for f_bm = 15% and f_bm = 25%. It can be noted that the model prediction agrees well with the clinical results. In addition to the standard once-weekly 70 mg dose, the BDG for a daily 10 mg dose is also represented, producing both treatments with very similar results.

FIGURE 8

FIGURE 8. Comparison of BDG vs. time plots of alendronate PK-PD1 model for (A) hip, that is, f_bm = 15% and (B) lumbar vertebra, that is, f_bm = 25% for a once-weekly 70 mg dose and daily 10 mg dose with the clinical results obtained by Rizzoli et al. (2002).

In Table 3, the results of the sensitivity analysis of the PK-PD1 model are presented. It can be seen that each constant affects at least one curve, having an important effect on the global error, except in the case where k_el,NCT is reduced by 10%, which does not change the fitting. The curves generated by the variations of each constant are included in the Supplementary Material.

TABLE 3

TABLE 3. Sensitivity analysis of the constants of the PK-PD1 model. The fitted constants of the PK-PD model were varied by ± 10%, one at a time, keeping the rest constant. $e_j^k$ is the error in each curve, where k = 1, 2 corresponds to f_bm = 15% and f_bm = 25%, respectively, and j = Chae, Kang, Khan, Riz indicates the curve. These error terms correspond to each of the summations divided by the corresponding mean in Eq. 30. E is the error as defined in Eq. 30.

If the influence of alendronate on the apoptosis and the resorbing capacity of osteoclasts is removed from the PK-PD model separately (models PK-PD2 and PK-PD3), the fitting errors are similar to those of the complete model (PK-PD1) (see Table 2). This indicates that both effects of alendronate have the same consequences, at least in the case simulated here with this model, as will be discussed later on.

If the saturation of the active subcompartment of BC is removed (model PK-PD4) or if no distinction is made between the active and inactive subcompartments, either with or without saturation (models PK-PD6 and PK-PD5, respectively), the fitting errors are higher (see Table 2) and therefore indicating that both model features, the distinction between subcompartments and the saturation, are important for the bone response. The plasma concentration curves are similar to those obtained with the model PK-PD1, but the fitting of the urine curves, particularly the long-term urine excretion, is significantly worse (see Figure 9). Regarding the bone gain, a difficulty is detected in adjusting the response for f_bm = 15% and f_bm = 25% simultaneously. In general, BDG is underestimated by models PK-PD4, PK-PD5, and PK-PD6 for f_bm = 15% and overestimated for f_bm = 25%, as can be seen in Figure 9.

FIGURE 9

FIGURE 9. PK-PD1, PK-PD4, PK-PD5, and PK-PD6 model results. (A) Alendronate short-term urinary excretion vs. time; (B) alendronate long-term urinary excretion vs. time; (C) BDG vs. time for hip, that is, f_bm = 15% and (D) BDG vs. time for lumbar vertebra, that is, f_bm = 25% for a once-weekly 70 mg dose.

Figure 10 shows the evolution of osteoclasts and osteoblasts population, bone resorption and formation rates, apparent density, and f_bm for model PK-PD1 in a 70 mg weekly treatment of alendronate and in the case of the hip. The results for the 10 mg daily treatment are very similar, and those of the lumbar spine show the same tendency. Alendronate treatment decreases both cell populations and thus reducing bone turnover rate. It also decreases the resorption capacity of osteoclasts k_res, and as a consequence, the formation rate is slightly higher than the resorption rate, thereby reversing the tendency of PMO. During the treatment, f_bm increases, but not enough to account for the BDG as can be deduced from the fact that the increment in apparent density is higher than that of f_bm, that is, the BDG is not only due to the reduction in porosity but also because of the mineralization of the tissue, produced by the decrease in bone turnover, similarly to what occurs in other antiresorptive treatments (Martínez-Reina and Pivonka, 2019).

FIGURE 10

FIGURE 10. PK-PD1 model results of 70 mg weekly treatment of alendronate in the case of the hip. Temporal evolution of (A) osteoblasts and osteoclasts concentration; (B) bone resorption and formation rates; (C) apparent density and f_bm. All the results are normalized with the values at the beginning of the treatment.

4 Discussion

The PK-PD model of alendronate proposed in this work was able to reproduce the short- and long-term pharmacokinetics reported in the literature (see Figure 7), as well as the BDG in a case of PMO for two different treatment protocols (10 mg administered daily and 70 mg administered weekly) for two different bone sites, that is, the lumbar spine and the hip (see Figure 8).

To model the pharmacokinetics of the drug, several models have been tested to conclude that the best option was a three-compartment model “in parallel.” Cremers et al. (2002) also used a three-compartment scheme in their pamindronate PK model (i.e., a different bisphosphonate); although in order to model its long-term treatment, the authors fitted the long-term urine excretion curves of the alendronate data obtained by Khan et al. (1997). These authors used an “in series” layout instead of our “in parallel” model. Differences may exist between the pharmacokinetics of pamindronate and alendronate, but, in any case, the fitting of both types of model (“in series” and “in parallel”) produced similar errors in our analysis, and we finally opted for the “in parallel” layout for two reasons: 1) the fitting error is slightly lower and 2) its biological basis is more solid. As Lin et al. (1991) experimentally demonstrated in animal models, alendronate is absorbed by both bone and non-calcified tissues, with a faster return to plasma from the latter compartment. Pillai et al. (2004) used a four-compartment model for ibandronate, although they stated that a three-compartment model was adequate. They introduced a fourth compartment on the basis to improve model performance. Other authors used models with a lower number of compartments (Sedghizadeh et al., 2013; Chae et al., 2014), but they did not fit the long-term experimental data. In our work, we conclude that one or two compartments are sufficient to predict the short-term behavior of bisphosphonates, but at least three compartments are required to also fit long-term experimental data.

The value obtained for bioavailability (see Table 2), F = 0.5%, is in agreement with the values found in the literature (Porras et al., 1999; Drake et al., 2008; Cummings et al., 2020).

The addition of the PD model did not significantly affect the results obtained with the PK model alone (see Figure 7). In the complete PK-PD model (PK-PD1), the return flow of alendronate from the BC to the CC is no longer independent of the osteoclast activity, that is, it is controlled by the elimination rate k_res Oc_a and therefore variable. However, the plasma concentration and the short- and long-term urine excretion have a similar goodness of fit compared to the purely PK model. For the first two quantities, the period of interest is significantly shorter than the time required to observe bone cellular changes during remodeling. For the long-term urine excretion, one would anticipate a significant influence. However, the osteoclast population is quickly stabilized and the same occurs with k_res (see Figure 10), therefore producing an approximately constant elimination rate, which is responsible for the urinary excretion rate in the long-term, slow but steady as seen in Figure 7.

The sensitivity analysis (see Table 3) showed that the variation of each constant has a significant effect on at least one of the curves. For the relatively small variation (10%), the PK constants do not have a significant effect on BDG, and the same occurs with the PD constants on the plasma concentration and urine excretion experimental curves. Another conclusion that can be drawn from the sensitivity analysis is the fact that the PK behavior (plasma and urine curves) is similar for f_bm = 15% and f_bm = 25% but not in the case of the PD behavior (BDG).

Our PK-PD model yielded very similar BDGs for the 2-year treatment period for both dose regimes (70 mg weekly and 10 mg daily), reproducing the clinical results obtained by Rizzoli et al. (2002) very closely in the lumbar spine and hip. This similarity between both protocols is also in agreement with the clinical results by Schnitzer et al. (2000). We identified that BDG can be due to two effects. First, to a decrease in bone turnover rate which affects osteoclasts more intensely. Their population falls less than that of osteoblasts (Figure 10C), but since the resorbing capacity (k_res) of the existing osteoclasts is also diminished by the drug, the resorption rate drops more acutely (Figure 10B). This fact makes f_bm rise slightly with the treatment. However, this rise cannot completely explain the predicted BDG, which was also influenced by the increased mineralization of the tissue produced by the decreased bone turnover, similarly to what occurs in other antiresorptive treatments (Martínez-Reina and Pivonka, 2019). The difference in BDG between sites is linked to the different loading and a dependence of the active alendronate on the bone-specific surface and f_bm. A noticeable decrease in the osteoclast population takes place during the first month of the treatment. Osteoblast population lags behind but also falls due to the coupling of bone cells in the first 6 months (see Figure 10). After this rapid decline, there seems to be a stabilization during the rest of the treatment. This behavior was experimentally reported by Rizzoli et al. (2002). Bone resorption markers (N-terminal telopeptide of type 1 collagen, NTx) and bone formation markers (bone-specific alkaline phosphatase, BSAP) exhibited strong decay in approximately the same times, faster for the NTx, and a subsequent stabilization, that is, however, more evident in the clinical results for the case of BSAP (Rizzoli et al., 2002).

As reported in the literature (Halasy-Nagy et al., 2001; Yu et al., 2018; Takagi et al., 2021), the role of alendronate in bone turnover was simulated in the proposed PK-PD model (PK-PD1) with a double mechanism of action: a weakening of the resorptive capacity of osteoclasts due to the disappearance of the clear zones and ruffled borders and an upregulation of their apoptosis. The influence of both mechanisms was evaluated by fitting two alternative models (PK-PD2 and PK-PD3) in which only one of them was considered alternatively. It was concluded that considering just one of those mechanisms (any of them) is enough to simulate the action of alendronate, as those alternative models produced similar errors in the fitting procedure (see Table 2). However, this result does not imply that any of them can be simply eliminated from the model, as it has been shown experimentally that both mechanisms take place. Probably, their consequences are similar to the clinical cases simulated here but might be different when longer periods of time are simulated or in more complex scenarios as discussed in the next paragraphs.

One simplification of BCPMs consists of assuming that bone turnover is continuous in time and space. The differentiation of mature osteoclasts and osteoblasts from their precursors occurs intermittently, and if there are no active osteoclasts at a given time point and bone site, then alendronate could not induce apoptosis at that site. On the other hand, only when the osteoclasts are being differentiated, alendronate could prevent the formation of clear zones and ruffled borders in developing osteoclasts. The two mechanisms are activated at different time scales, given that alendronate rapidly induces apoptosis (Takagi et al., 2021), but the differentiation of osteoclasts is a slower process. This would make them not totally equivalent in intermittent models.

Another limitation of the model is the use of a systemic PK model in conjunction with a local PD model. The clinical data of alendronate serum concentration and urine excretion used to fit the model are systemic. Indeed, the PK model depicted in Figure 4 is systemic itself. This means that the CC represents the serum of the whole body and the BC, the whole skeleton. However, the PD model was applied locally, at a certain RVE, in a specific bone site. Therefore, Eq. 15 and the last term of Eq. 21 imply a strong simplification of the model, which consists of assuming that the flux of alendronate from BC to CC at a given bone site is representative of what occurs at the skeleton on average. This concerns f_bm and Oc_a and implies that both variables are assumed to change globally as they do locally. A comprehensive PK-PD model should take into account the flux from every site of the BC following Eq. 15, with the specific values of f_bm and Oc_a at each point. Such a model would need to simulate the whole skeleton, but it must be disregarded for practical reasons.

If not the entire skeleton, at least one complete bone could be simulated with the proposed model. In these organ models, the concentration of alendronate could be heterogeneous with specific bone sites being saturated with the drug and others with low concentrations. The former would correspond to an active part of the BC, where the flux from CC to BC is hindered by the saturation, but the opposite flux would be favored if a resorption event occurred locally. On the contrary, the latter would correspond to the inactive part of the BC, where the drug is mainly accumulated and not released. Therefore, the distinction between active and inactive subcompartments would come naturally.

It must be noted that the predictions of models PK-PD4, in which the active subcompartment does not exhibit saturation, and PK-PD5, which ignores the distinction between the active and inactive subcompartments, almost coincide. This was expected, because the difference between both models is that the alendronate effect is proportional either to the active or total alendronate, respectively. The former is f_bm dependent, while the latter is not, but the variation of porosity is limited in the simulations carried out here. Thus, the active alendronate remains as a constant fraction of the total as long as saturation is not considered.

Figure 9 shows that the fit of the short-term urinary excretion and BDG curves, though not excessively bad, is worse than in the model PK-PD1, and the prediction of the long-term urinary excretion is rather poor. Therefore, in view of the results of models PK-PD4 and PK-PD5, which do not consider saturation, we can conclude that this feature is necessary. The need for saturation is supported by two further reasons. 1) In the clinical results of BDG (see Figure 8), a steep increase is observed at the beginning to slow down shortly thereafter, that is, the effect of alendronate is stronger during the first months, which reinforces the necessity of saturation. 2) According to Lin et al. (1991), alendronate would only stay in the blood for a very short time: some are quickly excreted via urine, and some are stored in the non-calcified tissues and the bone, but the return from the former to the blood is faster. So, we can deduce that alendronate is quickly deposited in the bone and returns to the blood via bone resorption, which occurs very slowly. According to this, the amount of alendronate accumulated within bone would increase constantly with each weekly dose, along with its effect on bone turnover which would also rise constantly, unless a saturation limit is reached. Clearly, saturation is needed to explain that the effect of the drug is reduced in the midterm.

Another question is how models PK-PD4 and PK-PD5 are able to achieve a reasonable fit for BDG without implementing saturation. The answer is that the mathematical fitting procedure is forced to completely readjust the pharmacokinetics, leading to some biological incoherences. First, as can be seen in Table 2, the BC absorption rate constant, k_BC, decreases by 4 or 5 orders of magnitude with respect to the standard model, which means that almost no alendronate enters the BC. Moreover, virtually none is stored because the elimination rate constant, k_el,BC, is higher. This prevents the accumulation of the drug within the bone matrix, and its effect on bone turnover is only achieved for enormous values of ${\mathrm{\Pi }}_{\mathit{rep}}^{\mathit{Ale}}$ and ${\mathrm{\Pi }}_{\mathit{act}}^{\mathit{Ale}}$ (see Table 2). Finally, although the fit of the urinary excretion curves is not very good with models PK-PD4 and PK-PD5, it would be even worse if the algorithm had not reduced k_el,NCT. Such reduction leads to higher disposal of the drug in the long term, though not from the BC but from the NCT, which would contradict the literature as the return to the serum is faster from the non-calcified tissues (Lin et al., 1991). In conclusion, the goodness of fit of these models is poor and worse, and they do not make any biological sense.

Once saturation proves necessary, the division of BC into an active and an inactive part also seems crucial, as is evident from the comparison between models PK-PD1 and PK-PD6. Both implement saturation but just of the active part in PK-PD1 and of the whole BC in PK-PD6. The latter is also unable to fit the clinical results, particularly the mid- and long-term urinary excretion because of the following. Saturation is needed to damp bone mass gain in the midterm, as commented before, but if saturation applies to the whole BC, then no more drug can enter the bone after saturation and must be excreted, thus failing to adjust the mid- and long-term excretion curves. In other words, the algorithm is not able to find a saturation limit that meets both requirements and either fails to predict BDG (if the saturation constraint is too loose) or fails to predict the excretion curve (if the saturation constraint is too tight). The latter is what can be seen in Figure 9 (top right, magenta line). The model is forced to reproduce the clinical results of BDG, but then the BC is saturated too soon and from that point on, all the administered drug is excreted at a constant rate. If only the active part of the BC is saturated (as in PK-PD1), BDG response may be damped while alendronate can continue to accumulate, now in the inactive part of the BC. This does not mean that this subcompartment cannot saturate, but that saturation is reached much later and probably only if the treatment is very long.

The PK-PD model of alendronate proposed in this work could provide a valuable tool to analyze the effect of alendronate in a large number of scenarios, particularly in PMO, as well as to design patient-specific treatments, including combinations of different drugs.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

JC-G, JM-R, PP, and RR-L conceived the study and planned the computational simulations. JC-G, JM-R, and RR-L developed the PK/PD model of alendronate. JC-G carried out the computational simulations. JC-G, JM-R, and PP wrote the manuscript. All authors discussed the results and reviewed and commented on the manuscript.

Funding

This study is part of the grant PID2019-106969R, funded by MCIN/AEI/10.13039/501100011033.

Acknowledgments

Dr. Pivonka gratefully acknowledges support from the Australian Research Council (ARC ITTC, IC190100020).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2022.940620/full#supplementary-material

Footnotes

¹It is important to note that if the PK model is not able to predict the concentration profile of the drug over time, the addition of a PD model is unlikely to improve the results.

²The rest goes to the inactive subcompartment, that is, Ale_BC,inact = Ale_BC − Ale_BC,act.

References

Acotto, C. G., Antonelli, C., Flynn, D., McDaid, D., and Roldan, E. (2012). Upper gastrointestinal tract transit times of tablet and drinkable solution formulations of alendronate: A bioequivalence and a quantitative, randomized study using video deglutition. Calcif. Tissue Int. 91, 325–334. doi:10.1007/s00223-012-9639-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashrafi, M., Ghalichi, F., Mirzakouchaki, B., and Doblaré, M. (2021). On the effect of antiresorptive drugs on the bone remodeling of the mandible after dental implantation: A mathematical model. Sci. Rep. 11, 2792. doi:10.1038/s41598-021-82502-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardoso, L., Fritton, S., Gailani, G., Benalla, M., and Cowin, S. (2013). Advances in assessment of bone porosity, permeability and interstitial fluid flow. J. Biomechanics 46, 253–265. doi:10.1016/j.jbiomech.2012.10.025

CrossRef Full Text | Google Scholar

Chae, J., Seo, J., Mahat, B., Yun, H., Baek, I., Lee, B., et al. (2014). A simple pharmacokinetic model of alendronate developed using plasma concentration and urine excretion data from healthy men. Drug Dev. Industrial Pharm. 40, 1325–1329. doi:10.3109/03639045.2013.819880

CrossRef Full Text | Google Scholar

Cocquyt, V., Kline, W., Gertz, B., Belle, S. V., Holland, S., DeSmet, M., et al. (1999). Pharmacokinetics of intrvenous alendronate. J. Clin. Pharmacol. 39, 385–393. doi:10.1177/00912709922007958

PubMed Abstract | CrossRef Full Text | Google Scholar

Cremers, S., Sparidans, R., Den, H. J., Hamdy, N., Vermeij, P., and Papapoulos, S. (2002). A pharmacokinetic and pharmacodynamic model for intravenous bisphosphonate (pamidronate) in osteoporosis. Eur. J. Clin. Pharmacol. 57, 883–890. doi:10.1007/s00228-001-0411-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Cummings, S., Santora, A., Black, D., Graham, R., and Russell, G. (2020). History of alendronate. Bone 137, 115411. doi:10.1016/j.bone.2020.115411

PubMed Abstract | CrossRef Full Text | Google Scholar

Drake, M., Clarke, B., and Khosla, S. (2008). Bisphosphonates: Mechanism of action and role in clinical practice. Mayo Clin. Proc. 83, 1032–1045. doi:10.4065/83.9.1032

PubMed Abstract | CrossRef Full Text | Google Scholar

Eastell, R., O’Neill, T., Hofbauer, L., Langdahl, B., Reid, I., Gold, D., et al. (2016). Postmenopausal osteoporosis. Nat. Rev. Dis. Prim. 2, 16069. doi:10.1038/nrdp.2016.69

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, J., Rodan, G., and Reszka, A. (2000). In vivo effects of bisphosphonates on the osteoclast mevalonate pathway. Endocrinology 141, 4793–4796. doi:10.1210/endo.141.12.7921

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleisch, H. (2000). Bisphosphonates in bone disease-from the laboratory to the patient. New York: Academic Press.

Google Scholar

Halasy-Nagy, J., Rodan, G., and Reszka, A. (2001). Inhibition of bone resorption by alendronate and risedronate does not require osteoclast apoptosis. Bone 29, 553–559. doi:10.1016/s8756-3282(01)00615-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernández, C., Beaupre, G., Marcus, R., and Carter, D. (2002). Long-term predictions of the therapeutic equivalence of daily and less than daily alendronate dosing. J. Bone Min. Res. 17, 1662–1666. doi:10.1359/jbmr.2002.17.9.1662

CrossRef Full Text | Google Scholar

Hernlund, E., Svedbom, A., Ivergard, M., Compston, J., Cooper, C., Stenmark, J., et al. (2013). Osteoporosis in the European Union: Medical management, epidemiology and economic burden. Arch. Osteoporos. 8, 136. doi:10.1007/s11657-013-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnell, O., and Kanis, J. (2006). An estimate of the worldwide prevalence and disability associated with osteoporotic fractures. Osteoporos. Int. 17, 1726–1733. doi:10.1007/s00198-006-0172-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, H., Hwang, S., Park, J., and Kim, C. (2006). HPLC method validation and pharmacokinetic study of alendronate sodium in human urine with fluorescence detection. J. Liq. Chromatogr. Relat. Technol. 29, 1589–1600. doi:10.1080/10826070600678308

CrossRef Full Text | Google Scholar

Khan, S., Kanis, J., Vasikaranand, S., Kline, W., Matuszewski, B., McCloskey, E., et al. (1997). Elimination and biochemical responses to intravenous alendronate in postmenopausal osteoporosis. J. Bone Min. Res. 12, 1700–1707. doi:10.1359/jbmr.1997.12.10.1700

CrossRef Full Text | Google Scholar

Lin, J., Chen, I. W., and Duggan, D. (1992). Effects of dose, sex, and age on the disposition of alendronate, a potent antiosteolytic bisphosphonate, in rats. Drug Metab. Dispos. 20, 473–478.

PubMed Abstract | Google Scholar

Lin, J., Duggan, D., Chen, I., and Ellsworth, R. (1991). Physiological disposition of alendronate, a potent anti-osteolytic bisphosphonate, in laboratory animals. Drug Metab. Dispos. 19, 926–932.

PubMed Abstract | Google Scholar

Martin, M., Sansalone, V., Cooper, D. M. L., Forwood, M., and Pivonka, P. (2019). Mechanobiological osteocyte feedback drives mechanostat regulation of bone in a multiscale computational model. Biomech. Model. Mechanobiol. 18, 1475–1496. doi:10.1007/s10237-019-01158-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, R. (1984). Porosity and specific surface of bone. Crit. Rev. Biomed. Eng. 10, 179–222.

PubMed Abstract | Google Scholar

Martínez-Reina, J., Calvo-Gallego, J., Martin, M., and Pivonka, P. (2022). Assessment of strategies for safe drug discontinuation and transition of denosumab treatment in PMO -Insights from a mechanistic PK/PD model of bone turnover. Front. Bioeng. Biotechnol. 10, 886579. doi:10.3389/fbioe.2022.886579

CrossRef Full Text | Google Scholar

Martínez-Reina, J., Calvo-Gallego, J., and Pivonka, P. (2021). Combined effects of exercise and denosumab treatment on local failure in post-menopausal osteoporosis–insights from bone remodelling simulations accounting for mineralisation and damage. Front. Bioeng. Biotechnol. 9, 635056. doi:10.3389/fbioe.2021.635056

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., García-Aznar, J., Domínguez, J., and Doblaré, M. (2008). On the role of bone damage in calcium homeostasis. J. Theor. Biol. 254, 704–712. doi:10.1016/j.jtbi.2008.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., and Pivonka, P. (2019). Effects of long-term treatment of denosumab on bone mineral density: Insights from an in-silico model of bone mineralization. Bone 125, 87–95. doi:10.1016/j.bone.2019.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Pillai, G., Gieschke, R., Goggin, T., Jacqmin, P., Schimmer, R., and Steimer, J. (2004). A semimechanistic and mechanistic population pk–pd model for biomarker response to ibandronate, a new bisphosphonate for the treatment of osteoporosis. Br. J. Clin. Pharmacol. 58, 618–631. doi:10.1111/j.1365-2125.2004.02224.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pivonka, P., Buenzli, P., Scheiner, S., Hellmich, C., and Dunstan, C. (2013). The influence of bone surface availability in bone remodelling — A mathematical model including coupled geometrical and biomechanical regulations of bone cells. Eng. Struct. 47, 134–147. doi:10.1016/j.engstruct.2012.09.006

CrossRef Full Text | Google Scholar

Porras, A. G., Holland, S. D., and Gertz, B. J. (1999). Pharmacokinetics of alendronate. Clin. Pharmacokinet. 36, 315–328. doi:10.2165/00003088-199936050-00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizzoli, R., Greenspan, S., Bone, G., Schnitzer, T., Watts, N., Adami, S., et al. (2002). Two-year results of once-weekly administration of alendronate 70 mg for the treatment of postmenopausal osteoporosis. J. Bone Min. Res. 17, 1988–1996. doi:10.1359/jbmr.2002.17.11.1988

CrossRef Full Text | Google Scholar

Sato, M., Grasser, W., Endo, N., Atkins, R., Simmons, H., Thompson, D., et al. (1991). Bisphosphonate action. alendronate localization in rat bone and effects on osteoclast ultrastructure. J. Clin. Invest. 88, 2095–2105. doi:10.1172/JCI115539

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnitzer, T., Bone, H., Crepaldi, G., Adami, S., McClung, M., Kiel, D., et al. (2000). Therapeutic equivalence of alendronate 70 mg once-weekly and alendronate 10 mg daily in the treatment of osteoporosis. Aging Clin. Exp. Res. 12, 1–12. doi:10.1007/BF03339822

CrossRef Full Text | Google Scholar

Sedghizadeh, P., Jones, A., LaVallee, C., Jelliffe, R., Le, A., Lee, P., et al. (2013). Population pharmacokinetic and pharmacodynamic modeling for assessing risk of bisphosphonate-related osteonecrosis of the jaw. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. 115, 224–232. doi:10.1016/j.oooo.2012.08.455

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, Y., Inoue, S., Fujikawa, K., Matsuki-Fukushima, M., Mayahara, M., Endo, Y., et al. (2021). Effect of nitrogen-containing bisphosphonates on osteoclasts and osteoclastogenesis: An ultrastructural study. Microsc. (Oxf). 70, 302–307. doi:10.1093/jmicro/dfaa073

CrossRef Full Text | Google Scholar

Ulrich, D., van Rietbergen, B., Laib, A., and Ruegsegger, P. (1999). The ability of three-dimensional structural indices to reflect mechanical aspects of trabecular bone. Bone 25, 55–60. doi:10.1016/s8756-3282(99)00098-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Valentín, J. (2002). Basic anatomical and physiological data for use in radiological protection: Reference values. Ann. ICRP 32, 1–277. doi:10.1016/S0146-6453(03)00002-2

CrossRef Full Text | Google Scholar

Yu, Z., Surface, L. E., Park, C. Y., Horlbeck, M., Wyant, G., Abu-Remaileh, M., et al. (2018). Identification of a transporter complex responsible for the cytosolic entry of nitrogen-containing bisphosphonates. eLIFE 7, e36620. doi:10.7554/elife.36620

PubMed Abstract | CrossRef Full Text | Google Scholar