INTRODUCTION
Bone remodeling is a physiologically regulated, dynamic process involving the continuous replacement of bone tissue and the adaptation of bone microarchitecture to changing mechanical and biochemical conditions [1, 2]. It is carried out by basic multicellular units (BMUs), composed of osteoblasts, osteoclasts, and osteocytes, whose activity is regulated by mechanical, hormonal, and molecular signals [2]. Remodeling allows for the maintenance of skeletal integrity, the repair of microdamage, and the adaptation to loads—both physiological and those induced by the presence of implant materials [3].
From the perspective of biomedical engineering, the ability to quantitatively model and predict spatial-temporal changes in bone density, structure, and composition is essential for optimizing the geometry of implants and endoprostheses [4], predicting the risk of osteolysis and post-implant resorption [5], simulating the effects of drugs on the skeletal system [6], and designing personalized therapies [1].
In recent decades, numerous numerical models have been developed to simulate the bone remodeling process. These are generally categorized into: phenomenological models, which describe remodeling as a material response to local mechanical stimuli (e.g. strain energy density, SED), while ignoring cellular activity and molecular pathways [4, 5, 7], and mechanobiological models, which explicitly implement biological regulatory mechanisms such as the RANK/RANKL/OPG pathway, osteocytic mechanotransduction, or the role of sclerostin and PTH [6, 8, 9, 10].
Despite growing biological knowledge, the translation of these insights into engineering practice remains limited. Phenomenological models, due to their low computational demands and ease of implementation (e.g. in FEM environments such as Abaqus, AN-SYS, FEBio), are widely used for analyzing implant integration [4, 5, 11, 12]. However, they neglect key biological phenomena, such as the impact of molecular mediators on osteoblast proliferation or variability in disease responses (e.g., osteoporosis, multiple myeloma) [13]. Mechanobiological models, while more biologically realistic, often remain conceptual—burdened by a high number of parameters, difficulty in calibration, and limited clinical validation [11, 14, 15, 16]
The objective of this article is a critical analysis of current numerical algorithms used in bone remodeling modeling, with particular emphasis on the range of biological processes captured, formal and computational complexity, engineering applications, translational barriers, and recommendations for future development.
BIOLOGICAL FOUNDATIONS OF BONE REMODELING - ASPECTS CAPTURED BY MATHEMATICAL MODELS GEOMETRICAL AND MATERIAL MODEL
The bone remodeling process is carried out by basic multicellular units (BMUs) composed of osteoblasts (bone-forming cells), osteoclasts (bone-resorbing cells), and osteocytes, which act as mechanosensors [2]. Mathematical models incorporate these cell types at varying levels of detail—from population-based variables (e.g., cell counts) to complex molecular interactions [8].
Osteocytes are considered the primary mechanosensory element. In most mechanobiological models, osteocytes initiate remodeling by releasing biochemical signals that modulate osteoblast and osteoclast activity [2, 9]. Models such as those by Pivonka [8] or Graham [26] implement this process as a biological activation function triggered by mechanical stimuli.
Osteoblasts form new bone tissue, and their population dynamics (e.g. proliferation, differentiation, apoptosis) are often described by equations influenced by factors such as TGF-β, sclerostin, or WNT/β-catenin signaling [6, 9]. Komarova [6] and Ji [17] are examples of such modeling approaches.
Osteoclasts resorb old bone tissue. Their formation and activity are typically governed by the RANK-RANKL-OPG signaling pathway. The RANKL/OPG ratio determines the resorption rate, a central mechanism in models by Komarova, Ayati, and Ryser [6, 13, 18].
The RANK/RANKL/OPG pathway is a core feature of nearly all mechanobiological remodeling models. RANKL, produced by osteoblasts and osteocytes, activates osteoclasts via the RANK receptor, while OPG acts as an inhibitor. Models describe the RANKL/OPG ratio as a control function for resorption rate [6, 8].
The WNT/β-catenin pathway regulates osteoblast proliferation and differentiation. It is activated by decreased levels of sclerostin—an inhibitor secreted by osteocytes. Graham and Pivonka incorporate this pathway to simulate osteoblast maturation [8, 9].
Sclerostin is a key inhibitor of bone formation, with expression levels regulated by local mechanical loading. In Martin’s model [19], a feedback loop is introduced linking tissue deformation to sclerostin levels, enabling dynamic simulation of biomechanical regulation.
PHENOMENOLOGICAL MODELS OF BONE REMODELING
Phenomenological models are a class of mathematical models in which bone remodeling is treated as a response of a deformable medium to local mechanical stimulus. They do not consider cellular mechanisms or biochemical regulation. These models are based on the theory of adaptive elasticity and define bone as a material that adapts its internal structure to mechanical conditions [5, 20].
Typically, these models use a governing equation of the form: dρ/dt = f(SED), where ρ is bone density and SED is strain energy density. Bone is considered to locally densify or resorb in response to mechanical demand. In numerical implementations, the function f is defined in a piecewise manner depending on whether the stimulus exceeds or falls below the physiological threshold [20].
Models developed by Huiskes [5], Beaupré [4], and Jacobs [12] laid the foundation for finite element-based remodeling simulations, used in predicting density changes around implants and prostheses [4, 5, 21]. These models introduce control functions that allow stable numerical solutions and facilitate implementation in FEM platforms (e.g. Abaqus).
More advanced phenomenological models incorporate time delay, asymmetry of formation/resorption rates, and non-linear control functions. However, they remain limited to mechanical stimulus as the driver of remodeling. Biological processes are not modeled, which limits applicability in scenarios involving pharmacological treatment or systemic diseases.
These models are successfully applied in predicting structural changes around endoprostheses, orthopedic implants, and in optimizing porosity in scaffolds. Their main advantages include low computational cost, implementation simplicity, and high numerical stability. Their main limitation is the lack of ability to simulate the effects of biological regulation, signaling, and cell-level interactions.
MECHANOBIOLOGICAL MODELS OF BONE REMODELING
Mechanobiological models represent the most advanced class of remodeling models. They combine mechanical stimuli (e.g., SED or strain) with biological regulatory processes, including cellular responses and signaling pathways. Rather than treating bone as an adapting material, these models simulate interactions among osteoblasts, osteoclasts, and osteocytes within the BMU [6, 8–10].
Foundational models by Lemaire and Komarova use coupled differential equations to describe populations of resorbing and forming cells [6, 22]. Central to these models is the RANK/RANKL/OPG signaling pathway, which controls osteoclast differentiation and survival. Models often incorporate effects of PTH, TGF-β, IGF-1, and interleukins [6, 8].
Graham’s model introduces osteocytes and the WNT/β-catenin pathway [9], while Martin and Scheiner link deformation to molecular signal expression [19, 23]. Pivonka’s comprehensive models integrate mechanical, biological, and temporal components [8]. Ayati and Ashrafi extend these models to include pathological and pharmacological influences [13, 24]. Kameo’s model simulates stress homeostasis mediated by osteocytic signaling [30].
These models are used to simulate drug action (e.g., denosumab), disease progression (e.g., osteoporosis, myeloma), and responses to molecular stimuli [6, 13, 24]. Despite their biological realism, they are used less frequently in engineering due to high complexity, numerous parameters, and difficulty of implementation in FEM environments. However, they offer promising potential for personalized medicine.
ADVANCED AND EMERGING MODELING APPROACHES
Recent advancements in bone remodeling theory have introduced several modeling strategies that go beyond classical phenomenological and mechanobiological paradigms. These advanced approaches aim to better replicate biological complexity by incorporating additional mechanisms such as microdamage accumulation, spatial diffusion of mechanical signals and multiphysics coupling involving fluid flow and biochemical transport.
Damage-based models represent an important class within this category. Unlike conventional strainor SED-driven formulations, these models posit that remodeling is regulated by the internal history of mechanical degradation. Addessi et al. [25] proposed a damage-dependent framework in which osteoclastic activation is triggered by the local accumulation of irreversible microdamage. This approach is particularly suitable for simulating long-term fatigue processes or pathological overload-induced resorption. Similarly, Dammak et al. [15] presented a computational scheme that couples evolution with adaptive bone turnover, enabling more accurate prediction of stress shielding and cortical thinning.
A different line of development focuses on the spatial nature of mechanotransduction. In diffusive-stimulus models, the mechanical signal responsible for initiating remodeling is assumed to spread through the bone matrix in a manner analogous to a diffusive field. Allena et al. [1] formulated such models using second-order partial differential equations that govern the transport of the remodeling stimulus across tissue regions. This approach mimics the biological reality of osteocyte network connectivity and fluid-based signaling, and has proven effective in capturing spatial heterogeneity in adaptation, especially in nonuniform anatomical sites such as the mandible.
Multiphysics models offer yet another level of biological fidelity by integrating mechanical deformation with interstitial fluid flow, ion diffusion, and solute transport. Cowin and Hegedus introduced a poroelastic theory of bone that links fluid pressure dynamics with mechanical stress redistribution, laying the foundation for modeling mechanotransduction in vascularized tissue [20]. This framework was further extended by Cowin and Weinbaum to account for solute-driven biochemical regulation, offering insight into remodeling in conditions of ischemia, osteoporosis, and implant-bone interface failure [16].
Giorgio et al. [26] developed an orthotropic continuum model with substructure evolution that interprets the primary mechanism behind Wolff’s law. Such models show considerable promise for simulating remodeling under multiaxial loading, pharmacological modulation, or systemic disorders affecting bone homeostasis.
Despite their conceptual richness, these advanced models share common challenges. Their parameter spaces are highdimensional and often difficult to calibrate, particularly due to limited availability of in vivo data. They also pose substantial computational costs, requiring solvers for coupled, nonlinear PDE systems. Nevertheless, their ability to replicate observed physiological and pathological behavior renders them powerful tools for hypothesis testing, implant optimization, and personalized simulation of bone adaptation.
To facilitate practical distinctions, Table 1 summarizes the principal model families by stimulus, biological fidelity, computational cost, and typical use.
Tab. 1.
Comparative summary of modeling families
| Model type | Main stimulus | Biological fidelity | Typical use/example |
|---|---|---|---|
| Phenomenological | SED/local strain | Low (no pathways) | Stress shielding; implant design [5] |
| Mechanobiological | Mechanical+ signaling | High (RANKL, WNT) | Drug/disease simulations [6, 8] |
| Damage-based | Accumulated microdamage | Medium | Fatigue-driven resorption [25] |
| Diffusive | Spatially diffused SED | Medium | Heterogeneous adaptation [1] |
| Multiphysics | Load+ fluid/solute transport | High–very high | Osseointegrati on [16, 20] |
In practice, the choice of model family entails distinct numerical burdens: phenomenological models offer unconditional stability and straightforward FE integration (e.g., density-update UMATs), whereas mechanobiological and multiphysics formulations require coupled ODE–PDE solvers, robust time-integration, preconditioned linear algebra, and careful parameter identifiability. For translational use, code availability, standardized datasets, and uncertainty quantification are as critical as mean accuracy.
APPLICATIONS AND COMPARISON OF NUMERICAL BONE REMODELING MODELS
Numerical models of bone remodeling are applied in implant design, prediction of structural changes, drug modeling, and metabolic disorder simulation. Model selection depends on the goal, available data, and computational resources.
Phenomenological models dominate in implant–bone interaction simulations due to their simplicity and compatibility with FEM. They are used in predicting bone density changes around hip and knee prostheses [4, 5, 21], dental implants [11,12, 26], and porous implant geometry optimization [27]. They do not require biological input data, which makes them accessible in clinical and engineering environments.
Mechanobiological models are used in analyzing drug effects (e.g. PTH, denosumab) [6], disease modeling (e.g. osteoporosis, cancer) [6, 13], and in molecular-level analysis of bone behavior [6, 17]. These models require complex calibration and are mostly used in research settings.
In practice, a dichotomy exists: phenomenological models dominate engineering applications, while mechanobiological models are more common in theoretical biology. Hybrid models that incorporate simplified biological mechanisms into mechanical frameworks may bridge this gap.
Patient-specific workflows for clinical translation. Image-based geometry (CT-derived surfaces/volumes) and regional density mapping enable subject-specific FE models, while calibration against individual follow-up data supports longitudinal prediction under changing loads or therapies. For clinical decision support, reporting parameter uncertainty, sensitivity, and robustness is as important as nominal accuracy, to ensure safe interpretation of model outputs in patient care.
Machine learning to complement mechanistic remodeling models. Data-driven methods can (i) provide surrogate models that accelerate FE simulations at design-space scale, (ii) enable Bayesian calibration/uncertainty quantification for parameters that are difficult to identify from sparse clinical data, and (iii) implement physics-informed learning (e.g., PINNs) to fuse governing equations with limited measurements. We note typical pitfalls—data shift, overfitting, and limited interpretability—and emphasize that ML serves to augment, not substitute, mechanobiological insight.
CALIBRATION AND VALIDATION IN MECHANOBIOLOGICAL MODELS
Calibration of mechanobiological models remains challenging due to parameter identifiability and the scarcity of prospective in vivo datasets. Representative strategies include (i) imaging-based mapping between CT/HU and elastic properties for patient-specific FEM, (ii) longitudinal follow-ups around dental or orthopedic implants comparing FE-predicted density/stress distributions with radiographic outcomes, and (iii) multiscale fits against observed adaptation patterns in vivo. Recent mandibular osseointegration studies illustrate imaging-based verification of FE predictions [11], while poroelastic and damage-diffusion formulations enable multiscale calibration of mechanotransduction and remodeling kinetics [16, 20, 25]. Establishing standardized validation protocols and data-sharing practices is essential to assess predictive utility for clinical decision support.
CONCLUSION
Numerical models of bone remodeling remain a key tool in analyzing skeletal processes, especially in biomedical, implantological, and pharmacological contexts.
Phenomenological models, due to low formal complexity, are widely used in tissue engineering and FEM simulations, although they do not capture complex biological mechanisms.
Mechanobiological models offer high biological realism, allowing simulation of drug and disease effects, but are challenging in terms of implementation and parameter calibration.
In engineering analyses, phenomenological models dominate, while mechanobiological models are used mainly in fundamental research.
The future of bone remodeling modeling lies in hybrid solutions integrating biological components into simplified numerical frameworks, along with greater integration of omics, histological, and imaging data.
In vivo validation studies are essential to assess the predictive efficacy of models in clinical contexts.
Recommendations for practice. For fast and robust structural predictions (e.g. stress shielding, early implant screening), phenomenological models are appropriate. For scenarios involving drug response, disease progression, or long-term remodeling, use mechanobiological or hybrid models with minimal-yet-salient biological regulators. For spatially heterogeneous adaptation or fatigue-related resorption, consider damage-based or diffusive-stimulus formulations. For vascular/metabolic coupling or osseointegration under ischemia, multiphysics models are preferred, provided that calibration and uncertainty reporting are feasible.
