In silico experiments of bone remodeling explore metabolic diseases and their drug treatment

Abstract

Bone structure and function are maintained by well-regulated bone metabolism and remodeling. Although the underlying molecular and cellular mechanisms are now being understood, physiological and pathological states of bone are still difficult to predict due to the complexity of intercellular signaling. We have now developed a novel in silico experimental platform, V-Bone, to integratively explore bone remodeling by linking complex microscopic molecular/cellular interactions to macroscopic tissue/organ adaptations. Mechano-biochemical couplings modeled in V-Bone relate bone adaptation to mechanical loading and reproduce metabolic bone diseases such as osteoporosis and osteopetrosis. V-Bone also enables in silico perturbation on a specific signaling molecule to observe bone metabolic dynamics over time. We also demonstrate that this platform provides a powerful way to predict in silico therapeutic effects of drugs against metabolic bone diseases. We anticipate that these in silico experiments will substantially accelerate research into bone metabolism and remodeling.

Fig. 1. In silico model of bone remodeling that incorporates mechano-biochemical couplings.

(A) Model of mechanosensing by osteocytes. Osteocytes produce mechanical signals S_ocy in response to a mechanical stimulus, defined as the modified equivalent stress σ_ocy (eq. S2 in Supplementary Methods S1.1), and transmit these signals to bone surface cells. S_r is a critical mechanical information that influences bone remodeling and is assumed to be the ratio of S_ocy to S_d, the latter being the average S_ocy over the region Ω. (B) Intercellular signaling for bone remodeling as incorporated into the bone remodeling platform (V-Bone). (C) Formulation of the spatial and temporal behavior of signaling molecules. The concentration of each signaling molecule ϕ_i is varied according to the reaction-diffusion equation, which includes production, degradation, diffusion, and reaction terms. (D) Probability of cell genesis, i.e., differentiation from precursor cells and proliferation and apoptosis for osteoclasts ($p_{\text{gen}}^{\text{ocl}}$, $p_{\text{apo}}^{\text{ocl}}$) and osteoblasts ($p_{\text{gen}}^{\text{obl}}$, $p_{\text{apo}}^{\text{obl}}$). These are regulated by the concentration of RANKL (RNL), Sema3A-Nrp1-PlxnA complex (SNP), sclerostin (SCL), and the mechanical information S_r, and can be described by Hill-type activator/repressor functions.

Fig. 2. In silico reproduction of bone adaptation to mechanical loading.

(A) Morphological changes by cooperative osteoclastic bone resorption (red) and osteoblastic bone formation (blue) in an inclined single trabecula (left) and a Y-shaped trabecula (right) under compressive load. Both trabeculae were compressed through elastic plates to attain 0.1% apparent strain along the z direction. (B) Three-dimensional model of a mouse distal femur reconstructed from microcomputed tomography images. This model was compressed to attain 0.1% apparent strain along the z direction, corresponding to the longitudinal direction of the femur. A cancellous bone cube with edge size 735 μm was selected as volumetric region of interest. (C) Morphological changes in trabeculae in the region of interest after 10 weeks of remodeling. A trabecula acquired the morphology suitable for supporting the load (red arrowhead), while a trabecula perpendicular to the loading direction was eroded (yellow arrowhead). (D) Measurement of the structural anisotropy of trabeculae in the region of interest using fabric ellipsoids based on the mean intercept length method. The lengths of the three principal semi-axes are denoted H_i, i = 1, 2, 3 (H₁ > H₂ > H₃). The degree of anisotropy, defined as H₁/H₃, increased from 1.28 to 1.43 after remodeling. For clarity, the fabric ellipsoid is displayed at twice its true size.

Fig. 3. In silico reproduction of osteoporosis and osteopetrosis caused by aberrant mechanical or biochemical conditions.

(A) Change in cancellous bone morphology after 5 weeks in a control model and an unloading model (in proximal view). In the unloading model, the applied uniaxial strain was ¹/₁₀ of that applied to the control model. Scale bar, 1 mm. (B) Enlarged views of cancellous bone in control and unloading models. Osteoclasts and osteoblasts on the trabecular surface are colored red and blue, respectively. Voxel size, 15 μm. (C) Quantification of changes in BV/TV, Oc.S/BS, and Ob.S/BS for 10 weeks in control (N = 5) and unloading models (N = 5). Oc.S/BS and Ob.S/BS are normalized by total bone surface. (D) Change in cancellous bone morphology for 10 weeks in an osteoporosis and osteopetrosis model (in proximal view). In these models, production of RANKL from the bone surface, exclusive of surface osteoclasts, was set to 1.3 and 0.7 times of that in the control model, respectively. Scale bar, 1 mm. (E) Quantification of changes in BV/TV, Oc.S/BS, and Ob.S/BS over 10 weeks in control (N = 5), osteoporosis (N = 5), and osteopetrosis models (N = 5).

Fig. 4. In silico perturbation of Sema3A to compare with corresponding in vivo experiments.

(A) Cancellous bone morphology in a mouse femur obtained by in vivo and in silico experiments on Sema3A-deficient mice. In the Sema3A-deficient model, production of Sema3A from the bone surface, exclusive of surface osteoclasts, was set to 0.5 times of that in the control model. Scale bar, 1 mm. (B) BV/TV and Tb.N as measured in vivo and in silico (N = 5). (C) Distribution of osteoclasts and osteoblasts on the trabecular surface immediately after starting simulation of control and Sema3A-deficient models. Voxel size, 15 μm. (D) Oc.S/BS and Ob.S/BS as measured in silico (N = 5). (E) Cancellous bone morphology in vivo and in silico in control and Sema3A-treated mice. Treatment with Sema3A was simulated by setting Sema3A production from the bone surface, exclusive of surface osteoclasts, to 1.5 times of that in the control model. Scale bar, 1 mm. (F) BV/TV and Tb.N as measured in vivo and in silico (N = 5). (G) Distribution of osteoclasts and osteoblasts on the trabecular surface after 5 weeks without treatment and immediately after starting Sema3A treatment. Voxel size, 15 μm. (H) Oc.S/BS and Ob.S/BS as measured in silico (N = 5). **P < 0.01; ***P < 0.005; NS, not significant, by Student’s t test.

Fig. 5. In silico prediction of the therapeutic effects of the osteoporosis drugs bisphosphonate (BP), anti-RANKL (RANKL-Ab), anti-sclerostin (SCL-Ab), and Sema3A.

(A) Cancellous bone morphology in a mouse femur modeled in silico without and with drug treatment. Upper panels show osteoporotic bones treated without and with drugs at high doses for 10 weeks. Lower panels are enlarged views. (B to D) Changes in (B) BV/TV, (C) Oc.S/BS, and (D) Ob.S/BS during drug treatment. (E) Rm.S/BS immediately after starting treatment with standard doses, and fraction of Oc.S/BS and Ob.S/BS in Rm.S/BS. (F) Apparent stiffness of cancellous bone along the loading direction after 10 weeks of drug treatment at standard dose. (G) Percentage changes in BV/TV and Oc.S/BS from the initial state when continuing or discontinuing anti-RANKL therapy. (H) Percentage changes in Ob.S/BS from the initial state when continuing bisphosphonate therapy or transitioning to anti-RANKL and anti-sclerostin therapy.