INTRODUCTION

Bone structure and function are maintained by homeostatic load-adaptive remodeling, which generates sophisticated bone microarchitecture to satisfy mechanical demands. This adaptive mechanism is the object of strong scientific and academic interest (1, 2). In addition, maintenance of load-bearing function throughout life is important to prevent bone fractures. Bone homeostasis can be disrupted by an imbalance between bone resorption and formation due to disuse or sex hormone aberrations, resulting in metabolic bone diseases such as osteoporosis (3, 4). Thus, it is indispensable to fully elucidate the underlying molecular and cellular mechanisms of bone metabolism and remodeling, from both scientific and clinical viewpoints.

Recent advances in molecular and cellular biology have helped identify multiple signaling pathways that regulate osteoclastic bone resorption and osteoblastic bone formation, as well as their relationship to mechanical stress (57). For example, genetic modification of signaling molecules in vivo has illuminated the molecular mechanisms of bone diseases (8, 9). These advances have also accelerated the development of molecularly targeted drugs against bone diseases (1012). However, the physiological or pathological status of bone as a system remains difficult to predict because of the interplay among bone cells and because of the complexity of relevant signaling networks.

To effectively prevent and treat bone diseases through a full understanding of bone remodeling regulated by mechano-biochemical couplings, computer simulation approachesthe so-called in silico approachesare of great significance. A large number of in silico researches on bone remodeling have been conducted by focusing on its mechanical aspect (13), and although they could reproduce adaptive changes of the bone microstructure to external loadings, the used in silico models were based on various phenomenological hypotheses regarding cellular mechanism. Increasing knowledge on cell-cell interaction via complex signaling pathways has motivated the development of in silico models that describe bone cell dynamics by explicitly taking into account the involved intercellular signaling (1416). These models allow theoretical evaluation of a biochemical aspect of bone remodeling. However, they cannot account for the relationship between spatially organized bone structure and the underlying cellular activities. Hence, a novel in silico model to investigate spatial and temporal behavior of bone remodeling that results from mechano-biochemical couplings is required.

We now enable simultaneous spatiotemporal observation of mechano-dependent intercellular signaling, bone cell dynamics, and bone morphological change through an in silico experimental platform (V-Bone) that mathematically models bone remodeling and links microscopic molecular/cellular interaction to macroscopic tissue/organ adaptation. The proposed in silico model was qualitatively verified from both mechanical and biochemical viewpoints by reproducing bone adaptation to mechanical loading and metabolic bone diseases. To quantitatively show the validity of the in silico model, in silico perturbation of a specific signaling molecule was conducted to compare with corresponding in vivo experiments. After quantitative validation, the in silico model was applied to predict the therapeutic effects of various drugs against osteoporosis. This platform is a revolutionary approach to fully and noninvasively explore bone remodeling dynamics over time, at scales ranging from the molecule/cell to the tissue/organ in a living body. The platform may also accelerate a paradigm shift in studies of bone metabolism and remodeling.

We propose an in silico model to investigate bone remodeling by incorporating mechano-biochemical couplings. Although bone remodeling is regulated by both local signaling factors and systemic hormones (17), to highlight osteocyte-driven bone remodeling as a local event, the in silico model is based on the assumption that mechanosensitive osteocytes buried in the bone matrix orchestrate osteoclastic bone resorption and osteoblastic bone formation via local intercellular signaling, without considering systemic hormonal changes. In addition, this model is based on the hypothesis that osteocytes regulate bone resorption and formation to achieve a locally uniform stress/strain state via bone remodeling (18, 19), which means that bone remodeling is susceptible to the local spatial variation of stress/strain in the bone tissue rather than their magnitude. Considering that osteocytes in the bone matrix are believed to be stimulated by interstitial fluid flow (20), which is driven by the gradient of fluid pressure instead of the fluid pressure itself, this hypothesis would be reasonable and is validated through a theoretical study (21).

In particular, mechanically stimulated osteocytes embedded in the bone matrix at position x produce the mechanical signal Socy (eq. S1), which is the product of the density of osteocytes ocy and the modified equivalent stress ocy, as shown in Fig. 1A (see Supplementary Methods S1.1). Through intercellular communication, the cell located on the bone surface at xsf integrates the local mechanical signals Socy within the neighboring region into Sd (eq. S3), which means the weighted average of Socy in . Ultimately, bone remodeling depends on the mechanical information Sr (eq. S5), a measure of local nonuniformity of stress defined by the ratio of Socy to Sd.

(A) Model of mechanosensing by osteocytes. Osteocytes produce mechanical signals Socy 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. Sr is a critical mechanical information that influences bone remodeling and is assumed to be the ratio of Socy to Sd, the latter being the average Socy 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 (pgenocl, papoocl) and osteoblasts (pgenobl, papoobl). These are regulated by the concentration of RANKL (RNL), Sema3A-Nrp1-PlxnA complex (SNP), sclerostin (SCL), and the mechanical information Sr, and can be described by Hill-type activator/repressor functions.

In response to mechanical stimuli, osteocytes activate or repress the activities of osteoclasts and osteoblasts through complex signaling cascades (see Supplementary Methods S1.2). An overview of intercellular signaling incorporated in the in silico model is presented in Fig. 1B. Sclerostin, a well-known mechanoresponsive protein in osteocytes that plays an important role in bone remodeling, inhibits osteoblastogenesis by binding to LRP5/6 and blocking canonical Wnt signaling and induces osteoblast apoptosis (6, 22). Production of sclerostin from osteocytes is reduced by mechanical loading (17, 22). On the other hand, the RANK/RANKL/OPG axis is primarily responsible for osteoclastogenesis. Osteoclast differentiation is induced by binding of receptor activator of nuclear factor- (RANK), which accumulates at the membrane of osteoclast progenitors, to RANK ligand (RANKL) produced by mesenchymal cells such as osteoblasts and osteocytes (6, 7, 23). In contrast, osteoprotegerin (OPG) released from mesenchymal cells inhibits osteoclast differentiation by sequestering RANKL (7, 12). Semaphorin 3A (Sema3A) inhibits osteoclast differentiation but promotes osteoblast differentiation by binding to a receptor complex consisting of neuropilin-1 (Nrp1) and one of the class A plexins (PlxnA) (24).

The spatial and temporal behavior of each signaling molecule is modeled as shown in Fig. 1C, where the concentration of signaling molecule i, i, varies according to the reaction-diffusion equation (eq. S8, see Supplementary Methods S1.3). The first, second, and third terms denote production, degradation, and diffusion of molecule i, respectively, while the last term describes the reaction of molecule i with molecule j, such as in ligand-receptor interaction (1416). We modeled mechano-biochemical coupling by describing the production rate of sclerostin PSCL as a monotonically decreasing function of the mechanical information Sr (eqs. S10 and S11), based on the experimental finding that Sost/sclerostin levels were reduced with increasing strain magnitude (25).

Bone remodeling is a cyclical process of bone resorption by osteoclasts and bone formation by osteoblasts (6, 7, 26). To express the initiation and termination of this cycle, the probability of cell genesis (i.e., differentiation from precursor cells and proliferation) pgeni and apoptosis papoi for bone surface cell i (i = ocl or obl) was modeled as a function of the concentration of signaling molecules (eqs. S24 to S27, see Supplementary Methods S1.4). As shown in Fig. 1D, the probability of osteoclastogenesis increases with the RANKL concentration but decreases with increasing Sema3A concentration. On the other hand, osteoblastogenesis increases with the Sema3A concentration but decreases with increasing sclerostin concentration. The probability of osteoblast apoptosis increases with the sclerostin concentration. An increase in mechanical information Sr was assumed to promote osteoclast apoptosis and inhibit osteoblast apoptosis.

By combining these in silico models (Supplementary Methods S1.1 to S1.4) with a voxel finite element method (FEM) for mechanical analysis (see Materials and Methods) (18, 27), we have constructed a unique and state-of-the-art in silico experimental platform (V-Bone) that incorporates mechano-biochemical coupling into bone remodeling.

Cancellous bone alters its trabecular orientation to coincide with principal stress trajectories, a phenomenon known as Wolffs law (2830). V-Bone enables observation of such mechanical adaptation in silico. To qualitatively verify the validity of the in silico model from a mechanical viewpoint, we reproduced bone adaptation to mechanical loading in a single trabecula and in cancellous bone spanning multiple trabeculae.

First, we simulated the adaptation of single trabeculae with two different configurations to compressive loading. A cylindrical trabecula with an inclined longitudinal axis was found to reorient parallel to the loading direction (Fig. 2A). In a Y-shaped trabecula, the branches moved toward each other. These results show functional adaptations in a single trabecula in response to external loads.

(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 Hi, i = 1, 2, 3 (H1 > H2 > H3). The degree of anisotropy, defined as H1/H3, increased from 1.28 to 1.43 after remodeling. For clarity, the fabric ellipsoid is displayed at twice its true size.

We then simulated the morphology of cancellous bone in a mouse distal femur subjected to physiological compressive loading using a model reconstructed from microcomputed tomography images (hereinafter called control model; Fig. 2B). Comprising multiple trabeculae within the inner cuboid region, most trabeculae acquired morphology suitable for supporting the load within 10 weeks (red arrowheads in Fig. 2C and movie S1). Several trabeculae perpendicular to the loading direction were also lost by bone resorption (yellow arrowheads in Fig. 2C). These results show that although individual trabeculae are networked in cancellous bone, they successfully adapt to imposed mechanical loads.

To quantify the adaptation in the region of interest, structural anisotropy was evaluated based on a fabric ellipsoid obtained by the mean intercept length method (18, 27). The direction of the three principal axes of the ellipsoid coincides with the principal directions of trabecular orientation, and their lengths indicate the characteristic lengths spanning bone and marrow space in the corresponding directions. Strikingly, the fabric ellipsoid stretched along the z direction as a result of 10-week remodeling (Fig. 2D), implying that cancellous bone acquired trabecular architecture totally parallel to the loading direction to satisfy the mechanical demand and suggesting functional adaptation at multiple trabeculae.

Together, the results indicate that by modeling complex intercellular signaling, V-Bone can reproduce bone adaptation to the mechanical loading, not only in a single trabecula but also in cancellous bone.

Osteoporosis, which is characterized by low bone mineral density and low bone quality, substantially reduces bone strength, leading to increased risk of bone fractures. The disease is triggered by low mechanical stress due to disuse (31) or by accumulation of factors that promote bone resorption, e.g., RANKL, due to sex hormone imbalance (57). On the other hand, osteopetrosis is one of inhered osteosclerotic disorders in which osteoclast dysregulation results in excess bone formation and bone hardening. Previously, we reported that conditional knockout of RANKL triggers osteopetrosis in mice (23). For qualitative verification of the in silico model from a biochemical viewpoint, we reproduced these metabolic bone diseases that include unloading-induced osteoporosis, as well as osteoporosis and osteopetrosis due to abnormal RANKL expression, by using multiple mouse femurs (N = 5).

We reproduced osteoporosis due to low mechanical stress, as observed in cases of extended bed rest and space flight (31). In particular, mouse femurs were simulated under low compressive loadings (hereinafter called unloading model) and compared with the control models (movies S2 to S5). In the unloading model, several trabeculae were lost around the central region of the femur (Fig. 3A), owing to excess bone resorption by osteoclasts at trabecular surfaces (Fig. 3B). Accordingly, the bone volume/tissue volume (BV/TV) ratio remarkably decreased in the first 2 weeks compared to that in the control model because of an increase in the osteoclast surface/bone surface (Oc.S/BS) ratio and a decrease in the osteoblast surface/bone surface (Ob.S/BS) ratio. Nevertheless, BV/TV plateaued after 2 weeks (Fig. 3C), indicating that cancellous bone adapts to the loss of external load within 2 weeks, at which point bone resorption and formation are again at equilibrium.

(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 1/10 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).

We also reproduced osteoporosis by up-regulating RANKL (hereinafter called osteoporosis model) (movies S6 and S7). In contrast to the unloading model, the osteoporosis model formed trabeculae throughout the femur (Fig. 3D, top). In addition, sustained activation of osteoclasts and slight inhibition of osteoblasts resulted in a gradual decrease in BV/TV over 10 weeks (Fig. 3E). These results imply that osteoporosis due to RANKL overexpression is characterized by chronic bone loss, while osteoporosis due to low mechanical stress is characterized by acute bone erosion (Fig. 3C). These observations are consistent with experimental data showing that BV/TV during bed rest or space flight decreases about 10 times faster than in primary osteoporosis (31).

Last, we reproduced an osteopetrotic state, which is characterized by abnormally high bone density, by down-regulating RANKL (hereinafter called osteopetrosis model) (movies S8 and S9). This model is characterized by increased trabecular thickness (Fig. 3D, bottom), with BV/TV monotonically increasing with time due to loss of RANKL-induced osteoclastogenesis (Fig. 3E).

Collectively, we have successfully simulated osteoporotic and osteopetrotic pathologies in silico, suggesting that V-Bone may reproduce a variety of metabolic bone diseases due to mechanical and biochemical determinants such as loss of mechanical stress and abnormal expression of signaling molecules.

Here, we describe an innovative approach to investigate the role of an essential signaling molecule in bone remodeling, in which the molecule of interest is perturbed in silico as is often done in vivo. Previously, mice deficient in Sema3A, a dual-function signaling molecule that inhibits bone resorption and promotes bone formation, were found to have a severe osteopenic phenotype due to osteoclast accumulation (24). Conversely, bone volume increases in mice treated with Sema3A, following loss of osteoclasts and accumulation of osteoblasts. We conducted in silico perturbation of Sema3A using multiple mouse femurs (N = 5) under the same conditions as in these in vivo experiments. Through quantitative comparison of the in vivo and in silico experimental results, the in silico model was validated.

Sema3A-deficient mice were modeled by down-regulating Sema3A (hereinafter called Sema3A-deficient model) and compared to the control model. Cancellous bone morphology in the Sema3A-deficient model was similar after 10 weeks of simulation to that obtained in vivo (Fig. 4A), with BV/TV and trabecular number (Tb.N) significantly smaller than those in the control model (Fig. 4B). In addition, the Sema3A-deficient model initially accumulated more osteoclasts at the trabecular surface to enhance bone resorption (Fig. 4, C and D). These results quantitatively resemble in vivo data (24).

(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 Students t test.

To investigate the therapeutic potential of Sema3A, bone remodeling was simulated for 5 weeks in the control model, followed by up-regulation of Sema3A for 5 weeks (hereinafter called Sema3A-treated model). The Sema3A-treated model generated thicker trabeculae than the control model after 10 weeks, as observed in vivo (Fig. 4E). The corresponding BV/TV and Tb.N values were also in close agreement with in vivo data (Fig. 4F). Immediately after Sema3A treatment, osteoblasts accumulated at the trabecular surface, as observed in vivo (Fig. 4, G and H).

Together, the data showed that in silico perturbation is a powerful way to clarify the effects of signaling molecules on bone dynamics at molecular/cellular and tissue/organ scales. Hence, such experiments may enhance the design of subsequent in vivo experiments and thus provide a novel approach to inspire and test new hypotheses regarding complex biological phenomena.

We propose a method to predict the therapeutic effects of various drugs against metabolic bone diseases in silico using V-Bone. We have now used this method to investigate the effects of dose, the resulting bone quality after drug treatment, and even the effects of different treatment regimens. In particular, we simulated the treatment of osteoporosis using bisphosphonate, anti-RANKL, anti-sclerostin, and Sema3A. Bisphosphonate, a current first-line therapy against osteoporosis, is specifically taken up by osteoclasts and is an inhibitor of bone resorption (11). Similarly, anti-RANKL potently inhibits bone resorption by suppressing osteoclastogenesis via RANKL (11, 12). Anti-sclerostin blocks binding of sclerostin to LRP5/6 and activates canonical Wnt signaling, thereby promoting bone formation and suppressing bone resorption (11, 22). Sema3A inhibits osteoclastic bone resorption and promotes osteoblastic bone formation (24). The effects of these drugs were modeled in V-Bone (see Supplementary Methods S1.5).

To predict the effects of patient-specific drug treatment, we simulated standard- and high-dose treatments (Fig. 5A) of one specific mouse femur, which are absolutely impossible to conduct in vivo. We assumed an idealized administration of each drug where the bioavailability is 100% and the plasma drug concentration is constant. In untreated osteoporotic models, BV/TV decreased from 18 to 9% after 10 weeks (Fig. 5B). At standard doses of all four drugs, BV/TV stabilized at about 15%. Standard doses also suppressed osteoclastogenesis (Fig. 5C). Whereas anti-sclerostin and Sema3A up-regulated osteoblastogenesis, bisphosphonate and anti-RANKL did not (Fig. 5D). These simulation results are consistent with the therapeutic effects reported in the in vivo experiments (32, 33). At high doses (threefold of the standard dose), antibodies to RANKL and sclerostin suppressed osteoclastogenesis (Fig. 5C), while anti-sclerostin and Sema3A enhanced osteoblastogenesis (Fig. 5D). Consequently, high doses of anti-sclerostin were the most effective in increasing BV/TV, while high doses of bisphosphonate exerted little influence on bone volume (Fig. 5B). Thus, therapeutic benefits gained from dose escalation substantially depend on the mechanism of action of the drug, highlighting the value of computational drug assessment in dose management.

(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.

In silico medication experiments enable analysis not only of bone quantity but also of bone quality, an important index for drug assessment. Although all four drugs stabilized BV/TV at almost the same level, the resulting bone quality varied, especially as assessed by repair of accumulated microdamage through remodeling (i.e., bone turnover rate) and by mechanical function to support external loads (i.e., bone mechanical integrity). Bone turnover rate was estimated as remodeling surface/bone surface (Rm.S/BS), also defined as the sum of Oc.S/BS and Ob.S/BS. Bone mechanical integrity was evaluated as the apparent stiffness of cancellous bone along the loading direction, a property that strongly depends on trabecular architecture (34, 35). Simulation results showed that administration of anti-sclerostin and Sema3A generates relatively high Rm.S/BS (Fig. 5E, left), mainly because of enhanced generation of osteoblasts (Fig. 5E, right). On the other hand, the apparent stiffness of cancellous bone after bisphosphonate therapy was lower than that after treatment with all other drugs (Fig. 5F). These results suggest that drugs that promote bone formation but inhibit bone resorption are more effective in improving both bone quantity and quality. The data also highlight that in silico experiments, unlike in vivo experiments, can simultaneously analyze cellular activities and mechanical properties for drug assessment.

Furthermore, in silico medication experiments provide a powerful way to predict the therapeutic effects of potential treatment regimens. For example, we simulated the following clinically relevant scenarios: discontinuation of anti-RANKL (36) and transition from bisphosphonate to anti-RANKL or anti-sclerostin (37). These scenarios were simulated to occur 5 weeks after treatment with the standard dose. Discontinuation of anti-RANKL decreased BV/TV at a constant rate (Fig. 5G, left) but rapidly increased Oc.S/BS, although the latter also gradually declined after peaking (Fig. 5G, right). These behaviors qualitatively coincide with clinical effects observed after discontinuation of anti-RANKL (36). Switching from bisphosphonate to anti-sclerostin increased Ob.S/BS to a larger extent than switching to anti-RANKL or retaining bisphosphonate (Fig. 5H). Together, the data suggest that V-Bone may potentially assist clinicians to devise previously untested treatment regimens before clinical trials.

We have developed a novel in silico experimental platform (V-Bone) to investigate spatial and temporal behavior of bone remodeling regulated by mechano-biochemical couplings, while previous in silico models of bone remodeling addressed bone structure/function and bone cell dynamics separately. The platform enables spatiotemporal observation and prediction of bone physiological and pathological conditions resulting from complex intercellular signaling. In conjunction with in vivo and in vitro experiments, in silico experiments provide a third avenue to explore bone metabolism and may thus accelerate research. Furthermore, we anticipate that V-Bone will prove valuable in clinical practice, such as in comprehensive drug assessment and formulation of effective treatment regimens.

The in silico model of bone remodeling was qualitatively verified from mechanical and biochemical viewpoints: We reproduced bone adaptation to mechanical loading (Fig. 2), as well as pathological bone states due to low mechanical stress and abnormal expression of signaling molecules (Fig. 3). To more rigorously validate the in silico model, we also demonstrated in silico perturbation of a specific signaling molecule, a standard in vivo technique in life science, and quantitatively compared the results with those from corresponding in vivo experiments (Fig. 4). In silico perturbation enables observation of the spatial and temporal dynamics of bone remodeling, which is difficult to achieve in vivo. Last, we applied the in silico model to predict the therapeutic effects of various drugs against osteoporosis and showed that in silico medication experiments provide a powerful way to assess the effects of drugs on bone cells and morphology in clinically relevant scenarios (Fig. 5). In all the in silico experiments conducted in the present study, mouse femurs were uniaxially compressed despite multiple loadings in vivo because of a lack of information on actual boundary conditions, which resulted in a unidirectional trabecular structure (Fig. 2, C and D). By considering more realistic loading conditions in the in silico model, which can produce various trabecular orientations (28), the reproducibility of the trabecular structure in response to mechanical loadings will be quantitatively corroborated by in vivo experimental data.

Measuring bone turnover markers and bone mineral density is the conventional noninvasive method to evaluate bone metabolic dynamics. Whereas this technique can measure temporal changes in the balance between bone resorption and formation, the resulting data do not include spatial information on bone morphology and cellular distribution. Although x-ray microcomputed tomography (38) can help bridge bone metabolism to the three-dimensional bone microstructure, live imaging of cellular behavior is difficult. Recently, intravital imaging of bone tissue has gained much attention as a new technique for real-time observation of spatiotemporal cellular activities (39), although it is only suitable for flat bone such as calvaria. In comparison to these experimental methods, V-Bone allows simultaneous spatiotemporal in silico observation and prediction of the distribution of signaling molecules, of bone cellular behaviors, and of bone microstructure.

An in silico experiment is an innovative way to explore molecular phenomena and thus will contribute invaluably to the progress of life science. The standard method to elucidate the role of a specific signaling molecule in a complex biological system is to test a research hypothesis in vivo, typically by perturbing the molecule of interest by techniques such as genetic manipulation. In contrast, we perturbed Sema3A in silico, a molecule that exhibits dual functions of inhibiting bone resorption and promoting bone formation, to emphasize the value of this approach. Bone morphometric data obtained in silico were quantitatively in good agreement with those obtained by corresponding in vivo experiments. These findings suggest that in silico perturbation may generate new research hypotheses that can then be tested in vivo, thereby accelerating hypothesis-test cycles to resolve outstanding research questions.

In silico medication experiments to predict the therapeutic efficacy of drugs against metabolic bone diseases are one of the promising clinical applications of V-Bone. Comprehensive in silico drug assessment at the preclinical phase of development will likely help clinicians determine the optimal strategy for drug administration and thus dramatically reduce the time and expense needed for large-scale clinical trials. In addition, in silico medication experiments will enable time-lapse evaluation of bone quality and bone quantity, especially because V-Bone uniquely predicts both cellular dynamics and tissue mechanical state in an individual patient. In the present in silico medication experiments, to focus on the relationship between mechanism of action of drugs and their therapeutic effects, we did not take into account the differences in bioavailability and biological half-life among the drugs. For clinical usage of V-Bone in the future, it is indispensable to incorporate these critical factors for pharmacokinetics that affect therapeutic efficacy. Thus, V-Bone may potentially enable personalized treatments for improving bone quantity and quality.

The concept of in silico experiments is greatly different from that of conventional computer simulations. In conventional computer simulations to capture the nature of complex phenomena through their replication, it has been considered that the number of parameters included in the in silico model should be kept as low as possible, and sensitivity analysis of these parameters can help us understand the essential characteristics of the phenomenon of interest. On the other hand, in silico experiments aim at observing the complex phenomena in silico as they occur in vivo to analyze the underlying mechanism and predict the events caused by arbitrary perturbation. Therefore, a model for in silico experiments is required to be constructed by taking into account the complexity inherent in the phenomena; hence, a large number of parameters are included in the in silico model (see tables S1.1 to S1.3). Some parameters that are difficult to determine directly from in vivo or in vitro experiments have to be set by heuristic methods. In addition, sensitivity analysis of all parameters included in an in silico model is almost impossible because of the huge degree of freedom. However, in the case of in silico experiments, parameter sensitivity analysis has the same meaning as in silico perturbation to investigate the effects of corresponding factors, which is a notable difference from conventional computer simulations. Once the proposed in silico model is validated through quantitative comparison with in vivo or in vitro experimental results, in silico perturbation of specific parameters in the in silico model, which is conventionally regarded as a parameter sensitivity analysis, holds promise for revealing an overlooked importance of unexpected factors.

Bone metabolism in our living bodies is regulated by many kinds of cells, such as hematopoietic stem cells and mesenchymal stem cells in the bone marrow, and related various signaling molecules. Furthermore, bone metabolism is coupled to a large biological system that includes endocrine, immune (40), and nervous systems (41). To highlight osteocyte-driven bone remodeling regulated by local signaling factors, in the present study, we explicitly modeled only osteoclasts, osteoblasts, and osteocytes, which are directly responsible for the change in bone volume, and several signaling molecules primarily relating to different functions. Despite these limitations, V-Bone was quantitatively validated through in silico perturbation of Sema3A (Fig. 4). This suggests that essential aspects of the actual complex molecular and cellular mechanism of bone remodeling can, to some extent, be represented by a reduced in silico model. Further expansion of V-Bone by incorporating other molecules or cells of interest will increase its prediction accuracy and expand the range of application. Thus, V-Bone is a promising framework, which potentially develops by including additional molecular and cellular mechanisms, to investigate the complexity inherent in bone remodeling and fully understand bone metabolism. We expect that experimental data about underlying molecular, cellular, and systems behaviors will accumulate exponentially in the near future. Accordingly, in silico experiments that integrate large data sets from in vivo and in vitro experiments will become more important as an alternative approach to investigate a wide range of molecular and cellular interactions. Incorporating quantitative in vivo and in vitro experimental data in the in silico models will enhance the validity of in silico experiments. We anticipate that V-Bone will accelerate bone metabolism and remodeling studies through comprehensive understanding of molecular, cellular, tissue, and organ dynamics.

Mechanical stress in bone tissue was analyzed by a voxel FEM (18, 27). Briefly, finite element models of mouse distal femurs (N = 5) were constructed from microcomputed tomography images. Each model was discretized using eight-node cubic finite elements with edge size 15 m. The bone was assumed to be homogeneous and isotropic, with Youngs modulus E = 20 GPa and Poissons ratio = 0.3. By using von Mises equivalent stress eq obtained through finite element analysis, the mechanical information Sr that regulates bone remodeling was determined (see Supplementary Methods S1.1).

The mechanical information Sr was coupled with intercellular signaling by influencing the production rate of sclerostin PSCL (eqs. S10 and S11, see Supplementary Methods S1.2 and S1.3). The reaction-diffusion equations (Fig. 1C and eq. S8) governing the spatial and temporal behavior of signaling molecules within the bone marrow were solved by an explicit finite difference method, in which the marrow space was discretized using the same voxel mesh as the FEM model. The effects of the drugs for osteoporosis were incorporated into the same equations (see Supplementary Methods S1.5).

The concentration of the specific signaling molecules and the mechanical information Sr determine the probabilities of cell genesis (i.e., differentiation from precursor cells and proliferation) pgeni and apoptosis papoi for osteoclasts and osteoblasts (i = ocl or obl) (Fig. 1D, eqs. S24 to S27, see Supplementary Methods S1.4). According to these probabilities, osteoclasts and osteoblasts are recruited on the bone surfaces or removed from them. The recruited osteoclasts and osteoblasts can alter the bone surface by resorbing old bone or forming new bone, respectively. To represent the changes in cancellous bone morphology, the level set method was used (42), because this method is capable of tracking the movement of individual trabecular surfaces (see Supplementary Methods S1.6). Cortical bone located near the outer surface of femurs was assumed to be not subject to morphological changes via remodeling.

Parameters used in in silico experiments are listed in tables S1 to S3. All the in silico experiments were carried out using an in-house code written in Fortran 90. The results were visualized with the open-source software ParaView (Kitware Inc.).

Acknowledgments: We thank Y. K. Kim, Y. Inoue, and K. O. Okeyo for discussions and comments. We also thank R. Terazawa for editing the figures and tables. Funding: This work was supported by Advanced Research and Development Programs for Medical Innovation (AMED-CREST), Elucidation of Mechanobiological Mechanisms and their Application to the Development of Innovative Medical Instruments and Technologies from Japan Agency for Medical Research and Development (AMED) (JP19gm0810003), the Acceleration Program for Intractable Diseases Research Utilizing Disease-Specific iPS cells from AMED (JP19bm0804006), Cross-Ministerial Strategic Innovation Promotion Program (SIP) (3D Design & Additive Manufacturing) from Japan Science and Technology Agency (JST), and VCAD System Research Program, RIKEN. Author contributions: Y.K. and T.A. designed the project and wrote the manuscript. Y.M. performed in silico experiments. M.H. and T.N. provided in vivo experimental data and collaborated with Y.K. and T.A. to plan the study and interpret the data. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data and source codes related to this paper may be requested from the authors.

See the article here:
In silico experiments of bone remodeling explore metabolic diseases and their drug treatment - Science Advances