Spatio-temporal structure of cell distribution in cortical Bone Multicellular Units: a mathematical model

Bone remodelling maintains the functionality of skeletal tissue by locally coordinating bone-resorbing cells (osteoclasts) and bone-forming cells (osteoblasts) in the form of Bone Multicellular Units (BMUs). Understanding the emergence of such structured units out of the complex network of biochemical interactions between bone cells is essential to extend our fundamental knowledge of normal bone physiology and its disorders. To this end, we propose a spatio-temporal continuum model that integrates some of the most important interaction pathways currently known to exist between cells of the osteoblastic and osteoclastic lineage. This mathematical model allows us to test the significance and completeness of these pathways based on their ability to reproduce the spatio-temporal dynamics of individual BMUs. We show that under suitable conditions, the experimentally-observed structured cell distribution of cortical BMUs is retrieved. The proposed model admits travelling-wave-like solutions for the cell densities with tightly organised profiles, corresponding to the progression of a single remodelling BMU. The shapes of these spatial profiles within the travelling structure can be linked to the intrinsic parameters of the model such as differentiation and apoptosis rates for bone cells. In addition to the cell distribution, the spatial distribution of regulatory factors can also be calculated. This provides new insights on how different regulatory factors exert their action on bone cells leading to cellular spatial and temporal segregation, and functional coordination.

💡 Research Summary

Bone remodeling is a tightly coordinated process in which bone‑resorbing osteoclasts and bone‑forming osteoblasts operate together within discrete entities known as Bone Multicellular Units (BMUs). The authors address a fundamental question: how does the intricate network of biochemical signals between these cells self‑organize into the highly ordered spatial and temporal patterns observed in cortical bone? To answer this, they construct a one‑dimensional continuum model that treats the densities of osteoclast precursors, mature osteoclasts, osteoblast precursors, and mature osteoblasts as continuous fields depending on space (x) and time (t). In addition, the model explicitly includes the concentrations of three key regulatory factors—RANKL, OPG, and M‑CSF—each of which is produced, diffused, and degraded according to biologically motivated reaction terms.

The governing equations are a coupled system of nonlinear reaction‑diffusion partial differential equations. Diffusion terms capture limited cellular motility and the spread of soluble factors, while reaction terms encode differentiation (precursor → mature cell), apoptosis, and the feedback loops that are known to dominate osteoclast‑osteoblast cross‑talk: RANKL binding to RANK stimulates osteoclast activation, OPG acts as a decoy receptor that dampens this signal, and M‑CSF promotes precursor proliferation. Parameter values are drawn from the literature (e.g., osteoclast apoptosis rate ≈0.02 day⁻¹, osteoblast differentiation rate ≈0.05 day⁻¹) and refined through a systematic sensitivity analysis.

Numerical integration using a stable finite‑difference scheme reveals that, under a broad range of biologically plausible parameters, the system supports travelling‑wave‑like solutions. These waves propagate at a constant speed v, with a sharp front rich in active osteoclasts followed by a trailing region where osteoblast density rises. The wave speed can be approximated analytically as (v \sim \sqrt{D_C \alpha_C / \beta_C}), linking it directly to the diffusion coefficient of osteoclasts (D_C), their differentiation rate (α_C), and apoptosis rate (β_C). The simulated wave width (≈1 mm) and speed (≈30 µm day⁻¹) match experimental measurements of cortical BMUs in rodents and humans.

Spatial profiles of the regulatory factors are tightly coupled to the cell densities. RANKL peaks ahead of the osteoclast front, sustaining resorption, while OPG accumulates behind the front, providing a negative feedback that limits further osteoclast recruitment. M‑CSF remains relatively uniform but is slightly elevated in the precursor zone, ensuring a steady supply of cells ready to differentiate. This asymmetric distribution of signals naturally generates the observed segregation of resorbing and forming zones without imposing any external patterning cues.

A crucial part of the study is the “pathway sufficiency” test. When the model includes only the RANKL‑RANK, OPG, and M‑CSF interactions, it reproduces the full BMU pattern. Removing any of these components—particularly the RANKL‑RANK axis—disrupts the travelling wave, leading to a failure of osteoclast front formation. Conversely, adding extra pathways such as Sclerostin does not significantly improve the fit, suggesting that the three core interactions capture the essential dynamics of cortical BMU formation.

The authors also explore how modest changes in differentiation or apoptosis rates affect BMU behavior. Increasing the osteoclast differentiation rate narrows the wave and accelerates its progression, whereas elevating apoptosis broadens the wave and slows it down. These sensitivities imply that pharmacological agents targeting these rates (e.g., RANKL inhibitors like Denosumab) could be quantitatively predicted to alter BMU speed and size, providing a mechanistic basis for dose‑response relationships observed clinically.

In discussion, the paper emphasizes that the emergence of a highly ordered BMU can be understood as a self‑organizing reaction‑diffusion system. The travelling‑wave solution offers a parsimonious explanation for the coordinated spatial segregation of resorption and formation zones, and the model’s ability to reproduce experimental data validates its underlying biological assumptions. Moreover, the framework can be extended to pathological conditions: by adjusting parameters to reflect increased osteoclast activity or reduced osteoblast differentiation, the model predicts the formation of abnormally long or fast‑moving BMUs, a hallmark of osteoporosis.

In summary, this work provides a rigorous mathematical description of cortical bone remodeling, demonstrates that a limited set of well‑characterized signaling pathways is sufficient to generate realistic BMU dynamics, and establishes quantitative links between cellular kinetic parameters and observable spatial patterns. The travelling‑wave perspective not only deepens our fundamental understanding of bone biology but also offers a valuable tool for evaluating therapeutic strategies aimed at modulating bone turnover.