Prediction of Fracture Risk From Early-Stage Bone Markers in Patients With Osteoporosis Treated With Once-Yearly Administered Zoledronic Acid.

The prevention of fractures is the ultimate goal of osteoporosis treatments. To achieve this objective, developing a method to predict fracture risk in the early stage of osteoporosis treatment would be clinically useful. This study aimed to develop a mathematical model quantifying the long-term fracture risk after 2 annual doses of 5 mg of once-yearly administered zoledronic acid or placebo based on the short-term measurement of bone turnover markers or bone mineral density (BMD). The data used in this analysis were obtained from a randomized, placebo-controlled, double-blind, 2-year study of zoledronic acid that included 656 patients with primary osteoporosis. Two-year individual bone resorption marker (tartrate-resistant acid phosphatase 5b [TRACP-5b]) and lumbar spine (L2-L4) BMD profiles were simulated using baseline values and short-term measurements (at 3 months for TRACP-5b and 6 months for BMD) according to the pharmacodynamic model. A new parametric time-to-event model was developed to describe the risk of clinical fractures. Fracture risk was estimated using TRACP-5b or BMD and the number of baseline vertebral fractures. As a result, the fracture risk during the 2 years was successfully predicted using TRACP-5b or BMD. The 90% prediction intervals well covered the observed fracture profiles in both models. Therefore, TRACP-5b or BMD is useful to predict the fracture risk of patients with osteoporosis, and TRACP-5b would be more useful because it is an earlier marker. Importantly, the developed model allows clinicians to inform patients of their predicted response at the initial stage of zoledronic acid treatment.

The primary goal of osteoporosis treatment is to prevent fractures. Bone mineral density (BMD), which accounts for 70% of bone strength, is considered the standard measure for the diagnosis of osteoporosis and the assessment of fracture risk. The effectiveness of therapies is evident long after the initiation of treatment because improvements in bone strength require time. In contrast, changes in bone turnover markers (BTMs) are observed within a few weeks after the initiation of treatments, including bisphosphonate. ,

Modeling and simulation consist of first constructing a mechanistic or semimechanistic mathematical model fitted to the observed data and then simulating several what‐if situations using the model to predict patient responses. This approach can be used to streamline drug development and personalize the drug dosage for each patient.

In the osteoporosis field, several previous studies have shown the importance and use of modeling and simulation. For example, Madrasi et al suggested that although both short‐term BTMs and long‐term BMD could predict fracture risk, BTMs may be more useful because they change faster than BMD. In addition, Bauer et al and Vasikaran et al reported that short‐term decreases in bone resorption markers could predict fracture risk. Ruixo et al developed a model to measure fracture risk, and Eudy‐Byrne et al related fracture risk to observed BMD data, but their models did not contain any drug‐specific information, such as doses or drug concentrations.

Annual intravenous administration of zoledronic acid, a third‐generation nitrogen‐containing bisphosphonate, has been approved in Japan for osteoporosis treatment. In a 2‐year randomized, placebo‐controlled study of 665 Japanese patients with primary osteoporosis, the once‐yearly administration of zoledronic acid (5 mg) reduced the risk of new morphometric vertebral and clinical fractures. Mori et al reported that BMD profiles up to 2 years could be simulated using patient background characteristics and the early response of tartrate‐resistant acid phosphatase 5b (TRACP‐5b). TRACP‐5b is one of the bone resorption biomarkers that negatively correlates with BMD and was demonstrated to associate with vertebral fracture in elderly women.

The aim of this study was to develop a mathematical model predicting long‐term clinical fracture risk after 2 annual administrations of 5 mg of zoledronic acid using the short‐term measurement of BTM or BMD in patients with osteoporosis.

Methods

Clinical Study Data

This analysis was performed on the basis of the secondary use of data obtained from a previous clinical study. The clinical study was conducted in compliance with the World Medical Association Declaration of Helsinki–Ethical Principles for Medical Research Involving Human Subjects and Good Clinical Practice. All data from this study were used in this analysis after anonymization.

Briefly, this was a multicenter, randomized, placebo‐controlled, double‐blind, parallel‐group, comparative, 2‐year study conducted in Japan with Good Clinical Practice standard. This trial enrolled 665 patients with primary osteoporosis who were randomly assigned to either the zoledronic acid 5 mg or placebo group in a 1:1 ratio. Randomized patients were administered the study drug once yearly by intravenous infusion (over 15 minutes). The 2‐year cumulative incidence of new morphometric vertebral fractures was 3.3% in the zoledronic acid group versus 9.7% in the placebo group (P = .0029). The numbers of any clinical fracture in the zoledronic acid and placebo groups were 24 (8.2%) and 52 (17.2%) for 24 months, respectively (P = .0014).

Among the full analysis set (FAS, n = 661), 5 patients were excluded due to a lack of TRACP‐5b data at 3 months, and thus the clinical data from 656 patients were used for this analysis. Although Nakamura et al used the incidence of new morphometric vertebral fractures as a primary end point, the occurrence of new morphometric vertebral fractures was too low and difficult to apply to develop a mathematical model. Thus, in this study, we analyzed the incidence of clinical fractures using our modeling approach based on a higher frequency. A clinical fracture was defined as a fracture confirmed on radiographs with clinically evident symptoms, such as pain on the back of the trunk or in the extremities, and thus did not include vertebral morphometric fracture without clinical symptoms. Clinical fractures were initially identified by physicians at each site and confirmed by radiographs or magnetic resonance imaging. The methods used to measure bone resorption markers, including TRACP‐5b, urinary N‐telopeptide of type I collagen (u‐NTX), and C‐telopeptide of type I collagen (CTX), and lumbar spine (L2‐L4) BMD were described in the previous report.

Model Development

Simulation of TRACP‐5b and BMD

TRACP‐5b level was measured at baseline before dosing and 3, 6, 12, 18, and 24 months after starting treatment. TRACP‐5b levels were simulated using the previously developed model. Briefly, the bone resorption markers were synthesized in a zero‐order manner with a constant Kin and eliminated with the first‐order rate constant Kout. The drug in the “effect site” was assumed to inhibit the synthesis of bone resorption markers (Figure 1). The model equations and their parameter estimates for the TRACP‐5b model are shown in Table S1. A 2‐year time course of TRACP‐5b for each patient was simulated with the individual parameter estimates calculated by an empirical Bayesian maximum a posteriori method using the first 2 observed data (baseline and 3 months after starting treatment).

Figure 1

Schematic description of the models linking TRACP‐5b and bone mineral density to fractures. A parametric time‐to‐event model linking bone resorption markers or BMD to fractures was used to quantify fracture risk. Effect site: hypothetical site where zoledronic acid is stored until the onset of the inhibitory effect. BMD, bone mineral density; KD, first‐order equilibrium rate constant; Kin, bone resorption marker production rate constant, estimated by Kout × individual baseline value of TRACP‐5b; Kout, bone resorption marker elimination rate constant; Ke0, equilibrium rate constant between TRACP‐5b and BMD; TRACP‐5b, tartrate‐resistant acid phosphatase 5b; Scale, ratio of change in BMD to change in bone resorption marker values.

BMD was measured at baseline and 6, 12, and 24 months after starting treatment. BMD levels were simulated using the same approach (Figure 1). The time profiles for BMD were assumed to be linearly affected by TRACP‐5b levels with the first‐order equilibrium constant Ke0, and they decreased with the same constant. The scale was used to adjust the change in bone resorption marker values according to the change in BMD. The model equation and its parameter estimates for the BMD model are shown in Table S1. A 2‐year time course of BMD for each patient was simulated by the empirical Bayes maximum a posteriori parameter estimates for each individual calculated using the first 2 observed data (baseline and 6 months after starting treatment).

Fracture Risk Model

Fracture risks were modeled by a parametric time‐to‐event model. The parametric hazard was assumed as exponential, Weibull, or log‐normal distribution, and the hazards were postulated to be related to the 2‐year TRACP‐5b or BMD profiles.

Covariate models were also tested using the following patient characteristics: sex, age, and baseline number of vertebrate fractures. Continuous covariates were modeled by the power model with standardization to their median values, and the power coefficients were estimated. Categorical covariates, such as sex and prior bisphosphonate use, were modeled in a relative effect manner.

The likelihood ratio test using a forward inclusion or backward elimination process was used to compare nested models and select the final model. The full model was developed by incorporating the statistically significant candidates after the forward inclusion process. Each candidate was then tested by removing them one by one to confirm their statistical significance. A minus twice log‐likelihood difference was used to assess the statistical significance of the candidate covariates. The P values for the forward inclusion process and backward elimination process were <.05 and <.01, respectively. The final model was generated using the remaining significant covariates. Akaike's information criterion (AIC) was also used to compare models, where a lower AIC indicates a better model.

The bootstrap method , , was used to estimate the standard errors for the estimates and evaluate the validity and robustness of the final model. Two hundred bootstrap replications were generated by random nonparametric resampling of the original data set with replacement. The final model was fitted repeatedly to the 200 data sets. All modeling and simulation in this study were performed using a nonlinear mixed‐effects model in Phoenix NLME 8.1 software (Certara LP, Princeton, New Jersey) with the Laplacian algorithm. Successful estimation was defined as the normal completion of the Phoenix software. The means of parameter estimates calculated from the successful estimations were compared with the final parameter estimates obtained from the original data set. A visual predictive check was also performed to visualize the model predictability. Two hundred random samples from the final model were simulated to give 90% prediction percentiles, which were compared to the observed Kaplan‐Meier curve.

Simulation Based on the Final Model

The effects of TRACP‐5b decreases from baseline to 3 months on the fracture risk were simulated using the developed model. Only patients with a baseline BMD T‐score <–2.5 (<–2.5 standard deviations below the average value for young healthy adults) were considered for this simulation because BMD T‐score <–2.5 is diagnosed as osteoporosis.

Results

This analysis used data from a total of 656 patients who had TRACP‐5b values at 3 months: 327 patients in the zoledronic acid group and 329 in the placebo group. Patient demographics and number of clinical fractures are shown in Table 1. No significant differences were seen for patient characteristics between the zoledronic acid and placebo groups.

Table 1

Patient Characteristics

	Zoledronic Acid, N = 327	Placebo, N = 329	Total, N = 656
Sex
Male, n (%)	21 (6.4)	19 (5.8)	40 (6.1)
Female, n (%)	306 (93.6)	310 (94.2)	616 (93.9)
Age, y	73.9 ± 5.3 (65‐88)	74.2 ± 5.4 (65‐87)	74.1 ± 5.4 (65‐88)
Weight, kg	52.3 ± 7.5 (35.7‐80.9)	52.1 ± 8.2 (32.1‐83.6)	52.2 ± 7.9 (32.1‐83.6)
Prior use of bisphosphonate
Never used, n (%)	295 (90.2)	302 (91.8)	597 (91.0)
Used with sufficient washout, n (%)	32 (9.8)	27 (8.2)	59 (9.0)
Baseline TRACP‐5b, mU/dL	416.6 ± 148.6 (143‐1240)	421.5 ± 159.7 (157‐1220)	419.1 ± 154.1 (143‐1240)
Baseline lumbar spine BMD (L2‐L4), g/cm2 a	0.680 ± 0.095 (0.36‐0.93)	0.674 ± 0.094 (0.39‐0.98)	0.677 ± 0.094 (0.36‐0.98)
Number of baseline vertebrate fractures (0/1/2/3/4/5/6)	29/165/86/32/11/3/1	35/160/84/37/12/1/0	64/325/170/69/23/4/1
Clinical fracture events, n (%)	24 (7.3)	52 (15.8)	76 (11.6)

BMD, bone mineral density; TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

Five patients were removed from this study due to a lack of TRACP‐5b measurements at 3 months. Values are expressed as means ± SD (range).

Lumbar spine BMD: n = 145 (zoledronic acid), 161 (placebo), and 306 (total).

Two‐year TRACP‐5b profiles for each patient taking zoledronic acid or placebo were simulated using the first 2 observed data at baseline and 3 months after dosing. The prediction showed good agreement with observation, although slight overprediction was seen after 18 months (Figure 2A and Figure S1).

Figure 2

Observed and simulated TRACP‐5b (A) or bone mineral density (B) values in patients receiving once‐a‐year zoledronic acid or placebo for 24 months. Blue symbols, observed data used for empirical Bayes individual estimations; red symbols, observed data not used in the parameter estimation. Lines indicate TRACP‐5b or BMD levels simulated by the empirical Bayes method using the model (Mori et al10). BMD, bone mineral density; TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

After confirming that a Weibull hazard model showed the lowest AIC value among parametric time‐to‐event models, the fracture risk model using Weibull hazard up to 2 years was developed with significant covariates of patients’ TRACP‐5b profiles as time‐varying covariate and number of prevalent vertebral fractures at baseline (Table 2). No other covariates, such as sex and age, were found to be significant. The simulated 90% prediction interval almost covered the observed Kaplan‐Meier fracture profiles, and the predictions were comparable to the observed fracture rate both for zoledronic acid and placebo groups (Figure 3). The final model showed that fracture risk increased when there was (1) a higher level of TRACP‐5b or (2) a larger number of baseline vertebral fractures.

Table 2

Parameter Estimates for the Time‐to‐Fracture Model Estimated by TRACP‐5b Profiles

Parameter	Estimate	CV%	2.5%CI	97.5%CI
γ	1.06	11.6	0.823	1.31
θ1	−9.37	9.10	−11.0	−7.69
θ2	0.781	29.2	0.333	1.23
θ3	0.351	28.9	0.152	0.550

CI, confidence interval; CV, coefficient of variation; TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

Weibull hazard model for fracture risk was used.

The number 400 was used for standardization as nearly equal to the median baseline value of TRACP‐5b (394.5).

NUMVFB indicates the number of baseline vertebral fractures.

θ2 > 0 indicates larger TRACP‐5b levels increase fracture risk.

θ3 > 0 indicates a larger number of baseline fractures also increases risk.

Figure 3

Observed and predicted fracture risk for 24 months. Simulated TRACP‐5b levels for each patient were used to calculate the fracture risk for 24 months. Red line, observed fracture risk; blue line and shade, predicted median with 90% prediction intervals. TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

An additional quantitative model to predict the fracture risk was developed using a subset of patients who had BMD observations (n = 306), thereby enabling us to compare 2 models using TRACP‐5b or BMD. The number of patients used for model development was smaller for BMD than for TRACP‐5b data (n = 656) because we used only patients with BMD measurements. The simulated BMD time course shown in Figure 2B indicates good agreement with observed BMD (shown in Figure S2) by using 2 data points—baseline and 6‐month values—into the developed model. The BMD model showed a slightly better fit with lower AIC (Table 3), where BMD effect was statistically significant and negative, indicating larger BMD improvement led to lower fracture risk. In addition, parameter estimates for the time‐to‐fracture model (BMD or TRACP‐5b model) using 306 patients are shown in Table 4 (BMD) and Table S2 (TRACP‐5b). The simulated fracture risk vs time is provided in Figure S3 and suggests good prediction coverage for the observed Kaplan‐Meier curve, although with relatively wide range because of the small number of fracture events.

Table 3

Comparison Among Models for Fracture Risk Quantification Using TRACP‐5b or BMD

Model	Predictor	N	AIC	Covariate Effect Estimate (95%CI)
1	TRACP‐5b	656	1442.318a	0.781 (0.333 to 1.23)
2	TRACP‐5b	306	657.425	0.521 (−0.158 to 1.23)
3	BMD	306	655.386	−0.115 (−0.250 to −0.018)

AIC, Akaike's information criterion (smaller is better); BMD, bone mineral density; CI, confidence interval; TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

Positive TRACP‐5b effect indicates larger TRACP‐5b levels increase fracture risk.

Negative BMD effect indicates larger BMD levels decrease fracture risk.

Model 1, which is the final model using TRACP‐5b and the data of 656 patients, is shown only for reference. The AIC for model 1 (1442.318) is not used for the comparison with other models because the number of patients affects the AIC value. Comparison should be made between models 2 and 3, which are based on the same number of patients.

Table 4

Parameter Estimates for the Time‐to‐Fracture Model Estimated by BMD Using the Data From 306 Patients

Parameter	Estimate	CV%	2.5%CI	97.5%CI
γ	1.06	13.3	0.851	1.41
θ1	−9.36	9.96	−11.6	−7.98
θ2	−0.115	50.4	−0.250	−0.018
θ3	0.352	58.5	−0.138	0.686

BMD, bone mineral density; CI, confidence interval; CV, coefficient of variation.

Weibull hazard model for fracture risk was used.

NUMVFB indicates the number of baseline vertebral fractures.

θ2 < 0 indicates a smaller increase in BMD also increases risk.

θ3 > 0 indicates a larger number of baseline fractures also increases risk.

Simulation Based on the Final Model

The effects of TRACP‐5b decreases (magnitude of 100, 200, or 300) from baseline to 3 months on the fracture risk were simulated using the model developed. Patients showing a large decrease in TRACP‐5b at 3 months after initiating treatment showed lower fracture risk during the subsequent 2 years (Figure 4).

Figure 4

Simulation of the effect of TRACP‐5b values on fracture risk. TRACP‐5b decrease from baseline to 3 months = 300 (red), 200 (green), and 100 (blue). Baseline TRACP‐5b = 400, number of baseline vertebral fracture = 1. Only patients with a baseline BMD T‐score < –2.5 (–2.5 standard deviations below the average value for young healthy adults) were included for this simulation. TRACP‐5b, tartrate‐resistant acid phosphatase 5b.

Discussion

To the best of our knowledge, this is the first study to develop a mathematical model to quantify fracture risks after the once‐yearly administration of zoledronic acid. The mathematical model developed here could accurately predict fracture risks using a short‐term (3‐month) measurement of the bone resorption marker TRACP‐5b or BMD at 6 months after the initiation of treatment. Although the BMD model resulted in a slightly better fit based on the AIC value, the TRACP‐5b model was considered more useful because it could predict future fracture risks at 3 months following treatment, whereas the BMD model required 6 months. It is highly desired to predict long‐term clinical responses in the early stage of drug treatments to help streamline decision making. This model allows clinicians to inform patients of their predicted response during the initial stage of zoledronic acid treatment. Therefore, the model using TRACP‐5b would be more useful.

Good agreement between predicted median with 90% prediction intervals and observed fracture frequency (Figure 3) showed that the predictions of fracture risk after zoledronic acid administration based on short‐term TRACP‐5b data using the developed model were valid. In this model, the values of 2‐year TRACP‐5b profiles, not the decrease from baseline, were used to quantify fracture risk. However, the simulation (Figure 4) showed that magnitudes in TRACP‐5b decrease were related to fracture risk for patients with a baseline BMD T‐score <–2.5, which was similar to our previous study. This suggests that the TRACP‐5b decreases at 3 months were important measurements to estimate the fracture risk of each patient for the next 2 years.

We used TRACP‐5b levels as a bone resorption marker to characterize fracture risk. We also tested the effect of either u‐NTX or CTX and reached a similar model (Table S3). TRACP‐5b is a specific marker for osteoclastic activity and the number of osteoclasts. In contrast, u‐NTX and CTX are markers of collagen breakdown. Thus, TRACP‐5b might have a more specific link to bone resorption, which may primarily reflect bone mineral content. Furthermore, the biological variability of TRACP‐5b is lower than that of u‐NTX and CTX, because TRACP‐5b is affected neither by renal dysfunction nor by food intake. Therefore, TRACP‐5b was considered the best predictor of fracture risk in our analysis.

The relevance of BTMs and fracture risk has not been fully established, and controversial results have been reported. Some studies showed relationships between decreases in short‐term bone resorption markers and the reduction in fracture risk. , , , , However, Bauer et al reported no significant relationship between short‐term changes in u‐NTX/creatinine or serum CTX and any fracture outcome. All of these controversial reports used descriptive statistical analyses or a simple linear regression model. To address this issue, a mathematical model with the ability to quantitatively link pharmacology to fractures would be a powerful tool. By developing the mathematical model described here, we have succeeded in applying short‐term BTM data to the prediction of long‐term fracture risk.

It has been reported that clinical vertebral fractures became more common in women. In this analysis, although sex was not found to be significant, the statistical power to detect sex difference might be low because male patients were very small relative to females (Table 1). The effect of sex on fracture risk prediction using TRACP‐5b should be clarified in a larger clinical study in the future.

To our knowledge, there are only a few studies that predicted fracture risk using statistical models. Ruixo et al predicted 3‐year vertebral fractures using a model for BMD time profiles. They constructed separate models for the active drug group and placebo group, and they did not directly analyze the drug effect, which is a clear contrast to our model. Eudy‐Byrne et al used observed BMD values (not a model prediction) to predict 10‐year clinical fractures. However, their model was also not a drug‐driven one.

In contrast, we have successfully developed a predictive model that can be used to predict the fracture risk of each patient. Therefore, the treatment effect of once‐yearly administered zoledronic acid can be predicted at an early stage (3 months) after the initiation of zoledronic acid treatment. We believe this prediction model will be valuable for the personalization of osteoporotic drug treatments.

When adherence to osteoporosis treatment is poor, as a result, the response of these patients to treatment is suboptimal. , , The once‐yearly administration of zoledronic acid is effective and is associated with minimal risk regarding low adherence, but a patient who takes zoledronic acid might not realize the effect in day‐to‐day life under a once‐a‐year medical treatment. It was previously reported that patients would be more likely to adhere to their medication if they received a positive message that they were responding to the therapy. , The prediction of fracture risk (Figure 4) could be used to inform patients of the effects of treatment and improve their satisfaction. An additional dosing of the once‐yearly preparation of zoledronic acid within a year is not approved when the predictive response might be suboptimal. Therefore, in case suboptimal response is predicted by our model, clinicians in charge are encouraged to monitor patients more intensively. And if undesirable symptoms such as bone loss might be observed, the clinicians can consider switching treatment from annual zoledronic acid to alternative drugs such as anabolic agents at an earlier stage of treatment.

A limitation of this analysis is the relatively short period (2 years) used to evaluate fractures. Although 2 years might appear short for evaluating osteoporosis treatment effects, there was a significant (P = .0014) decrease in any fracture event between zoledronic acid and placebo groups, suggesting that a 2‐year follow‐up was sufficient to evaluate fracture risks following zoledronic acid administration. Therefore, we believe that 2 years is not too short for the model‐based prediction of fracture risks.

Conclusions

We have developed a novel mathematical model to estimate future fracture risks after the administration of once‐yearly zoledronic acid. Fracture risk could be predicted based on TRACP‐5b or BMD profiles. During the early stage of treatment, 3‐month TRACP‐5b values are a highly useful marker to predict the fracture risk in 2 years. This model allows clinicians to inform patients of their predicted response at the initial stage of once‐yearly zoledronic acid treatment.

Conflicts of Interest

Y.M. and A.O. are employees of Asahi Kasei Pharma. Y.T. received grants from Asahi Kasei Pharma during the conduct of the study. H.K. and M.S. declare no conflicts of interest.

Author Contributions

Study design: H.K. and Y.T. Study conduct: H.K. and Y.T. Data collection: Y.M. and A.O. Data analysis: H.K. and YT. Data interpretation: H.K., A.O., M.S., and Y.T. Drafting article: H.K. and Y.T. Revising article content: H.K. and Y.T. Approving final version of article: H.K., Y.M., A.O., M.S., and Y.T. Y.T. takes responsibility for the integrity of the data analysis.

Additional supplemental information can be found by clicking the Supplements link in the PDF toolbar or the Supplemental Information section at the end of the web‐based version of this article.

Click here for additional data file.