Modeling Osteocyte Network Formation: Healthy and Cancerous Environments

Advanced cancers, such as prostate and breast cancers, commonly metastasize to bone. In the bone matrix, dendritic osteocytes form a spatial network allowing communication between osteocytes and the osteoblasts located on the bone surface. This communication network facilitates coordinated bone remodeling. In the presence of a cancerous microenvironment, the topology of this network changes. In those situations, osteocytes often appear to be either overdifferentiated (i.e., there are more dendrites than healthy bone) or underdeveloped (i.e., dendrites do not fully form). In addition to structural changes, histological sections from metastatic breast cancer xenografted mice show that number of osteocytes per unit area is different between healthy bone and cancerous bone. We present a stochastic agent-based model for bone formation incorporating osteoblasts and osteocytes that allows us to probe both network structure and density of osteocytes in bone. Our model both allows for the simulation of our spatial network model and analysis of mean-field equations in the form of integro-partial differential equations. We considered variations of our model to study specific physiological hypotheses related to osteoblast differentiation; for example predicting how changing biological parameters, such as rates of bone secretion, rates of cancer formation, and rates of osteoblast differentiation can allow for qualitatively different network topologies. We then used our model to explore how commonly applied therapies such as bisphosphonates (e.g., zoledronic acid) impact osteocyte network formation.

Introduction

Advanced cancers commonly metastasize to bone where they often disrupt the normal bone remodeling process (Roudier et al., 2003; Zhang et al., 2013). With the onset of various types of bone cancer, it is common for the bone remodeling process to be disrupted. Much previous work has been focused on macroscopic properties of the resulting bone, e.g., the osteoblastic (net bone formation) and osteolytic (net bone reduction) phenotypes. As bone is formed and resorbed cyclically, osteocyte networks can be morphologically malformed. A relatively unexplored area regarding cancerous bone formation is the study of osteoblast-to-osteocyte differentiation whilst concurrently taking into account network structure. Understanding the full nature of bone formation networks is now becoming extremely important, especially as osteocyte network structure is suspected to limit effectiveness of current anti-cancer therapy (Lerebours and Buenzli, 2016).

Myeloma and (benign) osteoma, osteocytes appear exceptionally spherical with shorter distorted dendrites that are reduced in number (Stinson, 1975; Eisenberger et al., 2008). An experimentally contrasting osteocyte network was observed with unregulated excessive cancerous growth in the presence of osteogenic sarcoma (Stinson, 1975). Broadly speaking, osteocytes within a cancerous microenvironment display either over or under developed phenotypes (see Table 1 and figures in Stinson, 1975) leading to “more connected” or “less connected” networks.

TABLE 1

Table 1. Cancer type impact on osteocyte topology.

Perturbations in osteocyte-network organization can impact both fluid flow and diffusion of metabolites and thereby affect mechanosenzation and signaling (Knothe Tate et al., 2000; Steck and Knothe Tate, 2005; Kerschnitzki et al., 2013; Lerebours and Buenzli, 2016). The exchange of signaling molecules through the lacuno-canalicular network relate to: skeletal unloading, and fatigue damage (Jilka et al., 2013). The range of signaling molecules that have been detected is vast and arise in the regulation of bone mineralization (Hesse et al., 2015), and many other organs¹. Studies have reported that high-density networks correlate positively with high bone quality (Kerschnitzki et al., 2013).

The functional role of these different network topologies is unclear. Experimental works only reveal a snapshot of the communication between the osteocytes within bone, and the osteoblasts on the bone surface. However, one can use mathematical modeling as a tool for investigation. Marotti et al. suggested that osteoblasts are incorporated into the osteocyte network by mature osteocytes extending their dendrites toward the osteoblast layer (Marotti et al., 1995; Marotti, 2000; Kamioka et al., 2001; Palumbo et al., 2004; Pazzaglia et al., 2010). Experimentally, sclerostin has been stained for and observed within the cancer structures between osteocytes and osteoblasts (Poole et al., 2005). These studies highlight the important roleosteocytes play in osteoblast function and differentiation.

Therapeutically, zoledronic acid is frequently used to treat metastatic breast cancer (BCA) where pathological bone is formed with lower densities of osteocytes per unit volume; zoledronate then helps recover the number of osteocytes by inhibiting osteoclasts—although it is not clear if the network structure is restored.

In this paper, we have developed a stochastic agent-based model to investigate how cancer cells regulate osteocyte behavior and bone formation. We consider osteocyte network formation building on an earlier model of osteocyte generation (Buenzli, 2015), which did not account for network structure. In said model, osteocytes are located within a growing domain that represents the presence of mineralized bone; osteoblasts are located on the surface of the bone substrate. There are two constitutive processes: (i) the bone surface moves with a speed proportional to the surface density of osteoblasts; and (ii) osteoblasts differentiate into osteocytes. We add to the model by allowing for extra structure relating to the osteocyte's canicular network. We allow for osteocytes to extend dendrites toward the osteoblast layer to signal osteoblast differentiation (we assume these processes occur concurrently). Osteoblasts are also allowed to proliferate and move along the bone surface.

With this model, we aim to link osteocyte density and network structure to biological quantities such as: the rate of osteoid secretion; the rate of osteocyte network formation; and the rate of preosteoblast proliferation. In particular, we investigate how stimulatory or inhibitory network dependent signals influence osteoblast-to-osteocyte differentiation and lead to different osteocyte network properties in newly formed bone.

1. Materials and Methods

1.1. Histological Analysis of Osteocyte Density

Specimens for analysis were derived from mice that were intratibially inoculated with either saline (Control) or breast cancer (BCa, PyMT cell line) as described previously under University of South Florida IACUC approved protocols (R2238 and R1762-CCL) (Araujo et al., 2014; Tauro et al., 2017). Mice were harvested for analysis prior to breach of the cortical bone (day 21). Separately, BCa-PyMT inoculated mice were treated with a bisphosphonate (zoledronate) over the course of the study period (1mg/Kg, subcutaneously thrice weekly) (Araujo et al., 2014; Tauro et al., 2017). Subsequent to tissue collection and isolation, bones were decalcified, processed and paraffin embedded. sections (5μm) were generated, rehydrated and then stained with either Gomori's Trichrome or Hematoxylin & Eosin using standard procedures.

We estimated the area of visible trabecular bone within a pathology image, and subsequently count the number of osteocytes within this region. For full details of this Algorithm, see Appendix A. Unfortunately, the samples were not at a resolution high enough to determine network structure, however it was possible to estimate number of osteocytes per unit area (#osteocytes/mm²). These results are presented as box plots in Figure 1.

FIGURE 1

Figure 1. (a) Tukey box plot showing results to preliminary image analysis of osteocyte number density. Means plotted in black (see text for details). The far left histogram (marked †) is representative of a collection of histology slides imaged at a lower resolution. (b) Comparison between osteocyte sizes: (left) control, (right) bone under breast cancer protocol. Image size corresponds to 32 × 21μm. The smaller box plots are the osteocyte densities for each mouse using 3–8 slides per mouse. Each sample of the osteocyte number density is weighted by the quantity of visible bone area within the sample when determining the data statistics summarized in the plots. The larger box plots show the combined data for mice undergoing the same protocol.

From Figure 1, breast cancer pathological bone have significant lower osteocyte number densities when compared to healthy bone. Macroscopically, breast cancer is often osteolytic and suppresses osteoblast proliferation and maturation. When applied as therapy, the zoledronate treatment allows for recovery of osteocyte number density.

1.2. Mathematical Model

1.2.1. Mathematical Description

Our model is a stochastic agent-based model for the bone formation phase only. The agents in our model are osteocytes and osteoblasts occupying nodes within a spatial network. The osteoblast-to-osteocyte differentiation pathway is subdivided into eight phenotypic stages: (i) preosteoblast; (ii) preosteoblastic osteoblast; (iii) osteoblast; (iv) osteoblastic osteocyte; (v) osteoid-osteocyte (Type II preosteocyte); (vi) Type III preosteocyte; (vii) young osteocyte; and (viii) old osteocyte (Franz-Odendaal et al., 2006). Additionally, the secretion of bone occurs as two steps: first osteoid is deposited as a collageneous scaffold, and then mineralization occurs to confer strength. Stages (iv)–(vi) are cells after the deposition front but before the mineralization front, surrounded by a non-mineralized osteoid matrix (i.e., there is scaffold around them). Stages (vii)–(viii) are cells whose volume has depleted (endoplasmic reticulum and Golgi apparatus reduction) and are in mineralized bone. The diagram in Figure 2 shows the bone-formation step. Here we are only interested in the structure of a mature osteocyte network [stages (vi)–(viii)], so we avoid modeling the full biology intricacy (e.g., cell sub-classifications, proteins etc) for simplicity. We model mobile preosteoblasts (disconnected from the osteocyte network) in front of the bone deposition front that proliferate; stationary mature osteoblasts that secrete osteoid and are connected to the osteocyte network; and osteocytes that form dendrites with mature osteoblasts.

FIGURE 2

Figure 2. Diagrammatic illustration of bone formation. (A) Biological process. Lighter shades of blue indicate more differentiated cells. The lighter shade of pink indicates the deposition front, and the darker shade of pink indicates the mineralization front. The left panel occurs earlier than the right panel. Dendritic osteocytes (light blue) have dendrites that extend toward the osteoblast layer (dark blue). The osteoblasts secrete bone matrix. Osteoblast cells marked with “A” are signaled by the osteocyte network to differentiate into osteocytes. Osteoblast cells marked with “B” do not differentiate and stay on the outer bone surface. Osteoblast cells marked with “C” arrive at the bone front after differentiating from precursor osteoblasts (preosteoblasts) (this figure is adapted from Franz-Odendaal et al., 2006). (B) Model representation of biological process. In a small time step, the following events can occur: (i) bone secretion; (ii) osteoblast differentiation; (iii) dendrite growth; (iv) osteoblast migration; and (v) osteoblast proliferation.

In our model, the spatial network has connections, between cells, representing the ability for two cells to communicate; physically this communication is mediated through dendrites.

The model consists of the following processes: (i) bone secretion; (ii) osteoblast differentiation; (iii) dendrite growth; (iv) preosteoblast migration; and (v) preosteoblast proliferation (see Figure 2). Agents of the same type (osteocytes or osteoblasts) follow the same rules, although each agent may have different properties, e.g., position, number of connections, in addition to type.

Bone secretion is carried out by osteoblasts secreting bone in a small region around themselves, orthogonal to the bone surface (in the normal direction). Osteoblasts can also become buried and differentiate into osteocytes. The rate of osteoblast differentiation may depend on if the osteoblast in question is in communication with the osteocyte network in bone.

For dendrite growth, osteocytes and osteoblasts create a communication channel (i.e., become connected in the network) at a rate that is a function of the distance between the two cells. Our model is consistent with the suggestions in Palumbo et al. (1990) that dendrites grow away from osteocytes within bone and toward the osteoblast layer on the bone surface.

Movement and cell division are included for pre-osteoblasts that are disconnected from the osteocyte network. Pre-osteoblasts move along the bone surface by means of a diffusive process; they are also able to proliferate and divide into create two daughter cells.

Our simulations are carried out in two dimensions, but they represent a slice from a three dimensional organ, in which the third dimension L_z is the typical distance between osteocyte centers (L_z ≈ 40μm), projected onto two dimensions (see Figure 7 in Appendix). The x-direction is the main direction of bone growth, and the volume occupied by the bone increases in time. We impose periodicity in the y-direction (orthogonal to bone growth) to avoid boundary effects.

In all the simulations in this paper, we consider bone formation after a cement line has been deposited. A cement line is a 1-5μm region of hypermineralized (and collagen deficient) bone (Skedros et al., 2005) is deposited after bone resorption to prepare surfaces for new bone formation. When this cement line is deposited, osteocytes from deep within the bone are not necessarily in communication with osteoblasts on the other side of the cement line (Qin et al., 1999).

Accordingly, we use the initial condition that there are no old osteocytes to communicate with (no initial network structure at the onset of bone formation) and there is an initial surface density of 6 × 10³mm⁻² pre-osteoblasts that do not have any network structure. We also specify that once an osteoblast is buried, a new osteoblast takes its place (so the total number of osteoblasts at any time is constant). This new osteoblast has no network structure. Therefore, one can interpret this configuration of the model as the scenario in which pre-osteoblasts are in abundance at the bone interface. New osteoblasts then move into the cell gaps on the bone surface as space becomes available (Figure 3).

FIGURE 3

Figure 3. Single simulation runs of bone growth with different parameter choices: (A) Healthy parameter set (as shown in Table 3); (B) Increased osteoid secretion (η → 2η); (C) Reduced dendrite growth (α → α/2); (D) Reduced osteoblast surface density ($\tilde{p}\to \tilde{p}/2$). Osteoblast are colored in blue and osteocytes (and their network connections) are in red; darker shades of red denote osteocytes that were buried earlier in time. All simulations are shown after 365 days.

1.2.2. Simulation and Analysis

Individual simulation runs of our stochastic agent-based model were performed with a fixed time step Monte Carlo algorithm (see Appendix B). However, a single realization is unrepresentative of the stochastic process, so many simulations must be carried out to work out statistics of the process along with parameterization. This can be computationally demanding, especially as the system grows in size.

As an alternative to large numbers of simulations, we can utilize mean field equations. We have shown previously (Taylor-King et al., 2017) that in the limit as the number of nodes in the network increases, one is able to derive a mean-field partial integro-differential equation for the expected number of nodes at a particular position in the domain connected to a fixed number of cells. By solving the mean-field equations, one can also calculate the degree distribution of the nodes of the network. The validity of the mean-field assumption is discussed in Taylor-King et al. (2017).

Simulations of the stochastic mathematical model suggest that the system is approximately homogeneous in y, the direction orthogonal to bone growth (see Figure 3). Under the assumption of homogeneity in y, the mean-field equations depend only on x, the direction of bone growth. The solutions of these equations provide spatial profiles of the densities of osteocytes and their network structure perpendicularly to the bone surface (see Appendix C).

Our mean-field equations take the form of two coupled hierarchies of integro-partial differential equations (see Equations 23–24 in Appendix). These equations can be solved to indicate how many particles one could expect to find within a small region of space are, connected to k other cells (referred to as degree k) at time t. We denote by v_k the surface density of osteoblasts (number per unit area, [v_k] = mm⁻²) connected to k osteocytes on the bone surface; and w_k the number density osteocytes (number per unit volume, [w_k] = mm⁻³) connected to k osteocytes or osteoblasts with position x at time t within the bone. Therefore, we write the total density of all osteoblasts (regardless of degree) as $p={\sum }_k{v}_k$ and total density of osteocytes (regardless of degree) as $q=\sum _k{w}_k$. The average number of osteocytes connected to an osteoblast is the average degree of osteoblasts, i.e.,

$\begin{array}{ll}{\langle k\rangle }_{\text{Ob}}=\frac{{\displaystyle {\sum }_{k=0}^{\infty }}{kv}_k}{p};& (1)\end{array}$

and the average number of cells (osteoblasts and osteocytes) in communication with an osteocyte is the average degree of osteocytes, i.e.,

$\begin{array}{ll}{\langle k\rangle }_{\text{Ot}}=\frac{{\displaystyle {\sum }_{k=0}^{\infty }}{kw}_k}{q}.& (2)\end{array}$

The full equations for v_k and w_k are derived in Appendix C and solved in Appendix D. In the present work, we only consider an uninterrupted bone formation process so that the density of osteocytes q at point x corresponds to that generated by the terminal differentiation of osteoblasts when the bone deposition front was at x.

The mean-field equations admit a traveling wave solution, corresponding to steady bone growth if the osteoblast proliferation rate is such that the surface density of osteoblasts is constant. Such a solution is useful both to explain the model and to investigate the qualitative effects of perturbations to parameters. In the following, we denote the traveling wave solution with a tilde.

2. Results

2.1. Selection of Differentiation Mechanism

We wish to investigate the effects of different model choices for the rate of osteoblast-to-osteocyte differentiation, D_k.

For a given surface density of osteoblasts, the volumetric density of inclusions embedded in a tissue during bone formation is determined by two dynamic processes: the rate of osteoblast terminal differentiation and the tissue growth rate (Buenzli, 2015). If the rate of osteoblast terminal differentiation is identical for all osteoblasts at a given location, i.e., D_k is independent of k ($Dk\equiv \widehat{D}$), the density of osteocytes generated is given by Buenzli (2015)

$\begin{array}{ll}\tilde{q}\equiv \frac{\widehat{D}\tilde{p}}{\tilde{\nu }}=\frac{\widehat{D}}{{\kappa }_{\text{form}}},& (3)\end{array}$

where $\widehat{D}$ is the terminal differentiation rate, i.e., the probability per unit time for a single osteoblast to become embedded as an osteocyte, $\tilde{\nu }={\kappa }_{\mathrm{\text{form}}}\tilde{p}$ is the normal velocity of the bone interface, and κform is the rate of bone deposition per osteoblast. If the rate of osteoblast terminal differentiation depends additionally on the number of connections k they have with osteocytes, the density of osteocytes generated sums up the contributions of all k-degree osteoblasts

$\begin{array}{ll}\tilde{q}\equiv {\displaystyle \sum _{k=0}^{\infty }}\frac{D_k{\tilde{v}}_k}{\tilde{\nu }}.& (4)\end{array}$

In contrast to Equation (3), the density of osteocytes given by Equation (4) depends explicitly on the population of osteoblasts, through the proportions $v_k/\tilde{p}$ of k-degree osteoblasts. Notice that this expression determines the density of osteocytes created at the moving bone deposition front. It does not account for processes that may subsequently affect osteocyte density such as osteocyte apoptosis, bone resorption and remodeling, which may remove and replace osteocytes subsequently (Lerebours and Buenzli, 2016).

To explore different differentiation mechanisms we alter both: the rate of dendritic growth, α, and the mechanism behind the rate of osteoblast differentiation, D_k.

2.1.1. Assuming No Network Influence: Degree-Independent-Rate Model

To account for the potential of network independent differentiation we assume that the number of osteocytes each osteoblast is in communication with does not impact the rate of osteoblast differentiation; we write

$\begin{array}{ll}{D}_k^{(\text{null})}\equiv \widehat{D},& (5)\end{array}$

for all values of k = 0, 1, …, ∞.

After parameterizing the system using experimental measurements found from literature (see Table 3), with this choice of osteoblast differentiation there remains two undefined (free) parameters in the model: the rate of dendrite growth $\widehat{\alpha }$, and the rate of osteoblast differentiation $\widehat{D}$; and we calibrate these parameters based on two further experimental observations: the mean degree of an osteoblast connectivity ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}=1$ (see Appendix E.6); and the number density of osteocytes $\tilde{q}=2.375\times 10^4{\mathrm{\text{mm}}}^{-3}$. By carrying out mathematical analysis in the traveling wave regime, we obtain formulas linking these quantities together (see Appendix D.2). We determine parameters as $\widehat{\alpha }=1.39\times 10^{-3}{\mathrm{\text{day}}}^{-1}$ and $\widehat{D}=2.59\times 10^{-3}{\mathrm{\text{day}}}^{-1}$.

If one changes the rate of osteoblast-to-osteocyte differentiation, $\widehat{D}$, we observe a linear relationship with the number density of osteocytes buried, $\tilde{q}$ (see Equation 3). Changing the rate of formation, $\widehat{\alpha }$, leads to a linear relationship with the mean osteoblast degree, ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}$, and the mean osteocyte degree, ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}$. These effects are decoupled — changing the rate of osteoblast-to-osteocyte differentiation, $\widehat{D}$, has no effect on the network structure, and changing the rate of dendrite formation, $\widehat{\alpha }$, has no effect on the osteocyte density, $\tilde{q}$. This will not be the case when osteoblast terminal differentiation is coordinated by the osteocyte network.

The assumption that osteoblast differentiation is independent of the osteocyte network is unlikely. First, as reviewed in Gohel et al. (1995), osteoblasts can be signaled by osteocytes to adhere to the mineral matrix and grow dendrites (subsequently differentiating) via the insulin-like growth factor 1 (IGF-1). Second, sclerostin secreted by osteocytes has been shown to act briefly as an inhibitory signal to prevent excessive osteoblast differentiation and allowing for coordinated osteocyte network formation (Poole et al., 2005). Note that in Poole et al. (2005), sclerostin was stained for and observable within the osteocyte's dendritic protrusions. These experiments, along with the three-dimensional scans taken by Kamioka et al. (2001), strongly suggest that osteocytes signal osteoblast differentiation through the extension of dendrites toward the osteoblast layer.

2.1.2. Modeling Network Effects

We now consider a range of models in which the pre-existing osteocyte network can either have an stimulatory or inhibitory effect on osteoblast differentiation. We consider differentiation rates of the form

$\begin{array}{ll}Dk=\{\begin{array}{ccc}\lambda & \text{if}& k=0,\\ \lambda +f(k)& \text{if}& k\ge 1,\end{array}& (6)\end{array}$

where λ is the network-independent rate of osteoblast differentiation and f = f(k) is the contribution to osteoblast differentiation for an osteoblast connected to k osteocytes. When f > 0 the network has an excitatory effect on osteoblast differentiation, and when f < 0 the network has an inhibitory effect on osteoblast differentiation. To prevent negative differentiation rates, we require that

$\begin{array}{ll}\lambda >0,\text{and}\lambda +f(k)\ge 0,\forall k\ge 1.& (7)\end{array}$

In the main text (see section 2.1.3 below), we consider only the case where f is a constant, f = γ. Other choices of f are considered in Appendix E, i.e., i.e., cumulative activation (f ∝ k) and diminishing activation (f ∝ 1/k). It should be noted that choices of f that have non-monotonic behavior (i.e., a local maximum/minimum exists) can lead to non-monotonic profiles in q.

2.1.3. Proposed Mechanism: Switch-Like Influence

Switch-like mechanisms are frequently found in biology (Cherry and Adler, 2000); at a cellular level this includes initiating mechanisms for proliferation and differentiation (Xia et al., 2006). For a switch-like osteoblast differentiation, we take

$\begin{array}{ll}{D}_k^{(\text{swt})}=\{\begin{array}{ccc}{\lambda }_{\text{swt}}& \text{if}& k=0,\\ {\lambda }_{\text{swt}}+{\gamma }_{\text{swt}}& \text{if}& k\ge 1.\end{array}& (8)\end{array}$

If the osteoblast on the bone surface is not in communication with any other osteocytes (k = 0) then it has an network-independent differentiation rate λ_swt. If osteoblasts are in contact with one or more osteocytes, then there is an induced differentiation rate λ_swt+γ_swt where γ_swt is the added contribution from the dendritic network.

Our technique for parameter identification with this family of terminal differentiation rates is based on the mean-field equations. We choose a network independent component of osteoblast differentiation such that ${\lambda }_{\mathrm{\text{swt}}}\ne \widehat{D}$, leaving 2 free parameters αswt and γ_swt. We then determine the free parameters by imposing that ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}=1$ and $\tilde{q}=2.375\times 10^4{\mathrm{\text{mm}}}^{-3}$ in the traveling wave regime. By making these assertions, if ${\lambda }_{\mathrm{\text{swt}}}<\widehat{D}$, the network contribution will always have an stimulatory effect on osteoblast differentiation; and if ${\lambda }_{\mathrm{\text{swt}}}>\widehat{D}$, then the network contribution will always inhibit osteoblast differentiation.

To investigate the role of intrinsic to extrinsic osteoblast differentiation rates, we assume that the network-independent differentiation rate λ_swt contributes half the total rate of osteoblast-to-osteocyte differentiation rate of the null model (so ${\lambda }_{\mathrm{\text{swt}}}=\widehat{D}/2$). We then determine that α_swt = 2.08 × 10⁻³day⁻¹ and γ_swt = 2.59 × 10⁻³day⁻¹. Another option would be for an inhibitory contribution, in which case we can set ${\lambda }_{\mathrm{\text{swt}}}=3\widehat{D}/2$, and then α_swt = 6.96 × 10⁻³day⁻¹ and γ_swt = −2.60 × 10⁻³day⁻¹.

2.1.4. Differentiation Mechanism Comparison

In Figure 4, we plot the osteocyte density profile, the mean osteoblast degree over time, and the mean osteocyte degree over time using Equations (23)–(25) in Appendix for 3 proposed choices of D_k: the null model, stimulatory network contributions, and inhibitory network contributions.

FIGURE 4

Figure 4. (Top) Osteocyte density profile [$q(t,x=\tilde{\nu }t)$], (Middle) mean osteoblast degree over time [${\langle k(x=\tilde{\nu }t)\rangle }_{\mathrm{\text{Ob}}}$], and (Bottom) mean osteocyte degree over time [${\langle k(x=\tilde{\nu }t)\rangle }_{\mathrm{\text{Ot}}}$] when solving Equations (23)–(25) in Appendix. The black line shows the null model, the blue line shows the stimulatory switch model, and the red line shows the inhibitory switch model.

Figure 4 shows that in all 3 models, it takes approximately 1 year (~0.7mm) to get to the steady-state desired osteocyte density. We see similar results for the mean osteoblast degree, except that inhibitory contributions to osteoblast differentiation require more time for the osteoblasts to reach the steady-state degree distribution. The main means for differentiation between the three models is the resulting mean osteocyte degree profile. The osteocyte density profile also differs at the onset, but slowly relaxes to the steady state profile. Keeping $\tilde{p}$, $\tilde{q}$, ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}$ fixed, stimulatory network contributions to osteoblast differentiation leads to more connected osteocyte networks with initially low numbers of osteocytes (higher final value of 〈k〉_Ot), and inhibitory network contributions lead to less connected networks with initially high numbers of osteocytes (lower final value of 〈k〉_Ot).

Capturing network structure and choosing a model that best reflects reality is difficult. For instance, the mean osteoblast degree profile takes a long period of time before settling to the steady-state mean osteoblast degree of ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}=1$. One assumption we have made throughout our model selection process was that the osteoblast surface density was approximately constant. It may be the case that this surface density changes over time (Eriksen et al., 1984; Parfitt et al., 1987). Were it the case that the surface density of osteoblasts was initially higher density before decreasing it would speed up the timescales until a steady-state traveling wave profile was reached.

2.2. Parameter Analysis

We explored a wide range of model parameters to further understand their impact on osteocyte network topology. For the switch-like model in section 2.1.3, we use the traveling wave analytic expressions described in Appendix D.2 to investigate how dependent variables $\tilde{q}$, $\tilde{\nu }$, ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}$, and ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}$ change due to a single perturbation in one of the independent variables η, $\tilde{p}$, κdiff, λ, γ and α, i.e., a sensitivity analysis, see Table 2. In addition to the dependent variables, we also include the number of dendrites per unit area $\tilde{{M}}=\tilde{q}{\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}/2$ which is a quantity that could be determined experimentally. Depending on whether osteocytes activate or inhibit the rate of osteoblast differentiation, we obtain different results for parameters λ, γ and α. For visualization purposes, we show a few stochastic simulation runs for a few specific examples in Figure 3; these examples use the switch-like model in described in section 2.1.3 with the inhibitory parameter configuration (${\lambda }_{\mathrm{\text{swt}}}=3\widehat{D}/2$, αswt = 6.96 × 10⁻³day⁻¹, γ_swt = −2.60 × 10⁻³day⁻¹).

TABLE 2

Table 2. Prediction summary.

Osteocyte density is mostly determined by the secretory rate and terminal differentiation rate, with some weak osteoblast density dependence via the frequency distribution of osteoblast degrees via equation (4). Table 2 confirms this as $\tilde{q}$ is dependent on the secretory rate (η) or the terminal differentiation rate (λ, γ, α), but not osteoblast density $\tilde{p}$. Contrasting these results to those given in Figure 4 and Table 2 shows intra-model variability by perturbing parameters, but Figure 4 shows inter-model differences by changing the mechanism for osteoblast differentiation.

Whilst many parameters change the density of osteocytes, only the rate of dendrite growth α and the rate of bone formation η can change the mean osteocyte degree in the steady state regime. A counter intuitive result also specifies that altering the rate(s) of osteoblast differentiation (parameters λ_swt and γ_swt) changes the osteocyte density (variable $\tilde{q}$) and the osteoblast degree (variable ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}}$), but not the degree of these newly formed osteocytes (variable ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}$). This is possible because osteoblasts are constantly being replaced in the steady state regime; therefore parameters λ_swt and γ_swt are changing the osteocyte density by modifying the mean time osteoblasts secrete bone before differentiation — without changing the mean osteocyte degree.

As we mentioned in 2.1.1, inhibitory network contributions to osteoblast differentiation is more likely in the presence of sclerostin, but stimulatory network contributions may also be possible via IGF-1. Determining which of these two mechanisms drives osteoblast differentiation will require further experimental work.

3. Discussion

We have presented a model for the formation of an osteocyte network and identified parameters for healthy bone formation. By perturbing parameters, one can investigate irregular bone formation and the resulting osteocyte topology changes. One can also then predict the driving differentiation markers that osteoblasts exhibit that would lead to these morphological changes. In the context of zoledronate therapy for breast cancer, we have used our model to propose how this commonly used clinical treatment impacts bone formation. For future experiments, we have suggested how measurable quantities link to underlying mechanisms.

The model proposed has some limitations, one could suggest many improvements to the model to improve our idealization of the osteocyte network; we discuss some of these below. Additionally, one might also want to consider the inclusion of chemical species representing proteins of interest, and the inclusion of osteoclasts to incorporate bone resorption. However, our model acts as a first step toward mathematically modeling osteocyte network formation, and avoids making overly specific assumptions on underlying mechanisms.

We now comment on how various aspects of our model compare to biological reality.

3.1. General Implications for Cancerous Bone Growth

A number of previous mathematical models have examined osteocyte density, but none of them have explored network structure. Graham et al. (2013), Moroz et al. (2006), and Wimpenny and Moroz (2007) give ordinary differential equation (non-spatial) models for cell populations; these include osteoblast, osteoclast, and osteocyte populations. Existing models of healthy bone remodeling (homeostasis) include spatiotemporal models (Ryser et al., 2009, 2010; Buenzli et al., 2012), but do not explicitly include osteocyte generation; these models have been adapted for the cancerous regime in Ryser et al. (2012). For TGFβ targeted therapy, modeling approaches have been used to optimize the treatment window of application (Cook et al., 2016). Mechanical focused methods capturing stresses and strains on the bone have also been explored (Rejniak and Anderson, 2011; van Oers et al., 2014). A general continuum modeling approach was also proposed in Buenzli (2015).

Our model-derived results show that osteocytes are either over differentiated (excessive dendrite growth) or underdeveloped (diminished dendrite growth). Additionally, we show that the osteocyte number density tends to decrease.

With experimental data that gives information on network structure (e.g., transmission electron microscopy, India ink histology stains), one should be able to approximately measure at least 2 of 3 quantities: the number of osteocytes present (quantitative estimate); whether the osteocytes are over-differentiated or underdeveloped (qualitative estimate); and finally the density of dendrites (quantitative estimate). Therefore, we should be able to compare a pathological bone slide to a (healthy) control slide and determine differences between: the osteocyte number density (q); the mean osteocyte degree (${\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}$); and the density of dendrites ($\tilde{{M}}$).

Furthermore, our model also makes testable predictions in regards to:

a.) If all 3 quantities (osteocyte number density, mean osteocyte degree, and dendrite area) have increased (resp. decreased), this corresponds to either: osteoblasts on the bone surface producing too little (resp. too much) osteoid when compared healthy bone; or that the rate of dendrite growth has increased (resp. decreased) in the excitatory switch model configuration.

b.) If the osteocyte number density has increased (resp. decreased) with an opposing decrease (resp. increase) in ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ot}}}$ or $\tilde{{M}}$, then this must correspond to the rate of dendrite growth changing but in the inhibitory switch model configuration.

c.) If the osteocyte density increases (or decreases), but the mean osteocyte connectivity remains constant, our model suggests this relates to a change the rate of osteoblast differentiation.

3.2. Implications for Zoledronate Treatment

Breast cancer is known to be osteolytic by promoting osteoclast mediated bone destruction. When applied as a therapy for breast cancer, the zoledronate treatment slows down bone resorption (along with other effects, Polascik and Mouraviev, 2008). This treatment is also associated with a recovery of osteocyte number density, see Figure 1.

Changing the proliferation ability of osteoblasts in the model changes the osteoblast surface density. This changes the quantity of bone produced per unit time, but not the network structure or steady state osteocyte density. In our model, changing the process of osteoblast/osteocyte maturation corresponds to changing either: the mechanism behind osteoblast differentiation (parameters λ, γ), or changing the rate of dendrite growth (parameter α).

In breast cancer (BCa), osteocytes have fewer dendritic connections to other osteocytes and osteocyte density is lower. Therefore the mean osteocyte degree of connectivity is reduced. To achieve a simultaneous decrease in osteocyte density and mean osteocyte degree in the model, one would either have to: decrease the rate of dendrite growth in an stimulatory switch model configuration; or increase the rate of bone secretion per osteoblast.

Given that the zoledronate treatment restores osteocyte density, we propose that future studies investigate how zoledronate acts on the osteoblast-to-osteocyte differentiation pathway. If bone treated with zoledronate has a healthy osteocyte network present, then one can assume zoledronate also targets (and restores) the rate of dendrite growth. However, if bone treated with zoledronate has a different network topology, one can conclude the differentiation mechanism is targeted, i.e., the network-independent or network-dependent rate of osteoblast differentiation has changed.

3.2.1. Triggering Osteoblast Differentiation

During a bone remodeling event, the total number of osteoblasts generated is far larger than the total number of new osteocytes generated (Parfitt, 1994). Pazzaglia et al. (2014) have estimated that only 1 in 67 osteoblasts become embedded in bone matrix as osteocytes over the depth of a single osteocyte. In our model using the steady state regime with the parameters from Table 3, over a length-scale of 5μm approximately 1 in 50 osteoblasts achieve terminal differentiation. The exact mechanisms behind this process are still poorly understood; they may involve physical processes such as burial by neighboring osteoblasts, or self-burial (Franz-Odendaal et al., 2006). It has also been suggested that subpopulations of osteoblasts are predestinated to become osteocytes (Marotti et al., 1992), and that this selection may be determined by the number of connections with osteocytes (Kamioka et al., 2001).

TABLE 3

Table 3. Model parameters, see Appendix E for more details.

One aspect ignored in our model is mechanotransduction. It is known that mechanical loads and fluid flow sheer stress can lead to greater dendrite growth (Zhang et al., 2006; Burra et al., 2010). After now posing our model, a pertinent question then remains as to the relative effect sizes between mechanical stimuli and microenvironmental signaling.

3.2.2. Osteocyte Degree Distribution

In the steady-state traveling wave regime, one can show that the node degrees of the osteocyte network are geometrically distributed when using either the null model (see section 2.1.1), or the switch-like proposed mechanism (see section 2.1.3). This effect comes from the difference equation structure shown in (23)–(25) in Appendix C.1.

In Kerschnitzki et al. (2013), a three-dimensional osteocyte network was studied, the topology of this network includes dendrites that do not connect to a second osteocyte, and connections that link between multiple osteocytes. Additionally, the nodes of their network included both osteocytes and the branching points of dendrites. Thus the functional communication network we studied is different from the lacuno-canalicular pore network and does not account for all types of communication redundancies that may exist. To incorporate all possible redundancies would require the inclusion of connections that do not connect to other nodes, and connections that exist between multiple nodes. We have some redundancy in our communication network in the form of multiconnections (multiple connections between two nodes). However in the limit of large networks, the probability of a multiconnection occurring in our model approaches zero. Even in the case of finite networks, this is very unlikely to occur in the model as ${\langle \tilde{k}\rangle }_{\mathrm{\text{Ob}}},\langle \tilde{k}\rangle \mathrm{\text{Ot}}={O}(1)$. Regardless of these model technicalities, it should be noted that the degree distribution of the three-dimensional scanned network in Kerschnitzki et al. (2013) was also shown to be geometrically distributed as derived from our model.

3.2.3. Orientation of Dendrites

It is clear from Figure 3 that we observe orientation of connections between older osteocytes and younger osteocytes or osteoblasts. This leads to an interesting question is as to whether one can configure our model to modify this orientation.

In brief, our model has allowed us to mechanistically explain osteocyte networks in the bone during bone remodeling and explain how these networks are impacted by cancer. Our approach identified, using experimental data, the parameters that characterize find new biology and derive new parameters A possible future direction in our work may be to modify our model to explore the functional difference between lateral connections, and connections perpendicular to the bone surface. As osteocytes have coordinated deposition, one may be able to explore whether coordinated terminal differentiation can occur as a function of lateral connections, burying a group of osteoblasts simultaneously.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

JT-K, PB, SC, and DB designed the research. JT-K produced the mathematical model under supervision by PB, SC, and DB. Biological analysis of bone sections was performed by JT-K and CL. The manuscript was written by all authors.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

JT-K acknowledges funding from the EPSRC (EP/G037280/1). PB gratefully acknowledges the Australian Research Council for a Discovery Early Career Research Fellowship (project number DE130101191). DB and CL acknowledge funding from the National Cancer Institute (U01 CA202958 and U01 CA244101). We thank Andrew Dhawan and Mason A. Porter for helpful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.00757/full#supplementary-material

Footnote

1. ^These include Receptor Activator of Nuclear Factor Kappa-B Ligand (RANKL), Vascular Endothelial Growth Factor (VEGF) (Jilka et al., 2013), Parathyroid Hormone (PTH) (Xiong et al., 2014), calcium ions (Ca²⁺) (Ishihara et al., 2012), and Sclerostin (Sapir-Koren and Livshits, 2014).

References

Araujo, A., Cook, L. M., Lynch, C. C., and Basanta, D. (2014). An integrated computational model of the bone microenvironment in bone-metastatic prostate cancer. Cancer Res. 74, 2391–2401. doi: 10.1158/0008-5472.CAN-13-2652

PubMed Abstract | CrossRef Full Text | Google Scholar

Buenzli, P. R. (2015). Osteocytes as a record of bone formation dynamics: a mathematical model of osteocyte generation in bone matrix. J. Theoret. Biol. 364, 418–427. doi: 10.1016/j.jtbi.2014.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Buenzli, P. R., Pivonka, P., and Smith, D. W. (2012). Bone refilling in cortical bone multicellular units: Insights into tetracycline double labelling from a computational model. arXiv preprint arXiv:1208.6075.

PubMed Abstract | Google Scholar

Buenzli, P. R., Pivonka, P., and Smith, D. W. (2014). Bone refilling in cortical basic multicellular units: insights into tetracycline double labelling from a computational model. Biomech. Model. Mechanobiol. 13, 185–203. doi: 10.1007/s10237-013-0495-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Buenzli, P. R., and Sims, N. A. (2015). Quantifying the osteocyte network in the human skeleton. Bone 75, 144–150. doi: 10.1016/j.bone.2015.02.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Burra, S., Nicolella, D. P., Francis, W. L., Freitas, C. J., Mueschke, N. J., Poole, K., et al. (2010). Dendritic processes of osteocytes are mechanotransducers that induce the opening of hemichannels. Proc. Natl. Acad. Sci. U.S.A. 107, 13648–13653. doi: 10.1073/pnas.1009382107

PubMed Abstract | CrossRef Full Text | Google Scholar

Cherry, J. L., and Adler, F. R. (2000). How to make a biological switch. J. Theoret. Boil. 203, 117–133. doi: 10.1006/jtbi.2000.1068

PubMed Abstract | CrossRef Full Text | Google Scholar

Cook, L. M., Araujo, A., Pow-Sang, J. M., Budzevich, M. M., Basanta, D., and Lynch, C. C. (2016). Predictive computational modeling to define effective treatment strategies for bone metastatic prostate cancer. Sci. Rep. 6:29384. doi: 10.1038/srep29384

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenberger, S., Ackermann, K., Voggenreiter, G., Sültmann, H., Kasperk, C., and Pyerin, W. (2008). Metastases and multiple myeloma generate distinct transcriptional footprints in osteocytes in vivo. J. Pathol. 214, 617–626. doi: 10.1002/path.2322

CrossRef Full Text | Google Scholar

Eriksen, E., Gundersen, H., Melsen, F., and Mosekilde, L. (1984). Reconstruction of the formative site in iliac trabecular bone in 20 normal individuals employing a kinetic model for matrix and mineral apposition. Metab. Bone Dis. Relat. Res. 5, 243–252. doi: 10.1016/0221-8747(84)90066-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Franz-Odendaal, T. A., Hall, B. K., and Witten, P. E. (2006). Buried alive: How osteoblasts become osteocytes. Dev. Dyn. 235, 176–190. doi: 10.1002/dvdy.20603

PubMed Abstract | CrossRef Full Text | Google Scholar

Gohel, A. R., Hand, A. R., and Gronowicz, G. A. (1995). Immunogold localization of beta 1-integrin in bone: effect of glucocorticoids and insulin-like growth factor i on integrins and osteocyte formation. J. Histochem. Cytochem. 43, 1085–1096. doi: 10.1177/43.11.7560891

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, J. M., Ayati, B. P., Holstein, S. A., and Martin, J. A. (2013). The role of osteocytes in targeted bone remodeling: a mathematical model. PLoS ONE 8:e63884. doi: 10.1371/journal.pone.0063884

PubMed Abstract | CrossRef Full Text | Google Scholar

Hannah, K., Thomas, C., Clement, J., Carlo, F. D., and Peele, A. (2010). Bimodal distribution of osteocyte lacunar size in the human femoral cortex as revealed by micro-CT. Bone 47, 866–871. doi: 10.1016/j.bone.2010.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Hesse, B., Varga, P., Langer, M., Pacureanu, A., Schrof, S., Männicke, N., et al. (2015). Canalicular network morphology is the major determinant of the spatial distribution of mass density in human bone tissue: evidence by means of synchrotron radiation phase-contrast nano-CT. J. Bone Miner. Res. 30, 346–356. doi: 10.1002/jbmr.2324

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishihara, Y., Sugawara, Y., Kamioka, H., Kawanabe, N., Kurosaka, H., Naruse, K., et al. (2012). In situ imaging of the autonomous intracellular Ca2⁺ oscillations of osteoblasts and osteocytes in bone. Bone 50, 842–852. doi: 10.1016/j.bone.2012.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Jilka, R. L., Noble, B., and Weinstein, R. S. (2013). Osteocyte apoptosis. Bone 54, 264–271. doi: 10.1016/j.bone.2012.11.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamioka, H., Honjo, T., and Takano-Yamamoto, T. (2001). A three-dimensional distribution of osteocyte processes revealed by the combination of confocal laser scanning microscopy and differential interference contrast microscopy. Bone 28, 145–149. doi: 10.1016/S8756-3282(00)00421-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerschnitzki, M., Kollmannsberger, P., Burghammer, M., Duda, G. N., Weinkamer, R., Wagermaier, W., et al. (2013). Architecture of the osteocyte network correlates with bone material quality. J. Bone Miner. Res. 28, 1837–1845. doi: 10.1002/jbmr.1927

PubMed Abstract | CrossRef Full Text | Google Scholar

Knothe Tate, M. L., Steck, R., Forwood, M. R., and Niederer, P. (2000). In vivo demonstration of load-induced fluid flow in the rat tibia and its potential implications for processes associated with functional adaptation. J. Exp. Biol. 203, 2737–2745.

PubMed Abstract | Google Scholar

Kristinsson, S. Y., Minter, A. R., Korde, N., Tan, E., and Landgren, O. (2011). Bone disease in multiple myeloma and precursor disease: novel diagnostic approaches and implications on clinical management. Expert Rev. Mol. Diagnost. 11, 593–603. doi: 10.1586/erm.11.44

PubMed Abstract | CrossRef Full Text | Google Scholar

Lerebours, C., and Buenzli, P. (2016). Towards a cell-based mechanostat theory of bone: the need to account for osteocyte desensitisation and osteocyte replacement. J. Biomech. 49, 2600–2606. doi: 10.1016/j.jbiomech.2016.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Logothetis, C. J., and Lin, S.-H. (2005). Osteoblasts in prostate cancer metastasis to bone. Nat. Rev. Cancer 5, 21–28. doi: 10.1038/nrc1528

PubMed Abstract | CrossRef Full Text | Google Scholar

Marotti, G. (2000). The osteocyte as a wiring transmission system. J. Musculoskelet. Neuronal Interact. 1, 133–136.

PubMed Abstract | Google Scholar

Marotti, G., Ferretti, M., Muglia, M., Palumbo, C., and Palazzini, S. (1992). A quantitative evaluation of osteoblast-osteocyte relationships on growing endosteal surface of rabbit tibiae. Bone 13, 363–368. doi: 10.1016/8756-3282(92)90452-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Marotti, G., Ferretti, M., Remaggi, F., and Palumbo, C. (1995). Quantitative evaluation on osteocyte canalicular density in human secondary osteons. Bone 16, 125–128. doi: 10.1016/8756-3282(95)80022-I

PubMed Abstract | CrossRef Full Text | Google Scholar

Moroz, A., Crane, M. C., Smith, G., and Wimpenny, D. I. (2006). Phenomenological model of bone remodeling cycle containing osteocyte regulation loop. Biosystems 84, 183–190. doi: 10.1016/j.biosystems.2005.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Palumbo, C., Ferretti, M., and Marotti, G. (2004). Osteocyte dendrogenesis in static and dynamic bone formation: an ultrastructural study. Anatom. Rec. A Discov. Mol. Cell. Evol. Biol. 278, 474–480. doi: 10.1002/ar.a.20032

PubMed Abstract | CrossRef Full Text | Google Scholar

Palumbo, C., Palazzini, S., Zaffe, D., and Marotti, G. (1990). Osteocyte differentiation in the tibia of newborn rabbit: an ultrastructural study of the formation of cytoplasmic processes. Cells Tissues Organs 137, 350–358. doi: 10.1159/000146907

PubMed Abstract | CrossRef Full Text | Google Scholar

Parfitt, A. (1994). Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone. J. Cell. Biochem. 55, 273–286. doi: 10.1002/jcb.240550303

PubMed Abstract | CrossRef Full Text | Google Scholar

Parfitt, A. M., Drezner, M. K., Glorieux, F. H., Kanis, J. A., Malluche, H., Meunier, P. J., et al. (1987). Bone histomorphometry: standardization of nomenclature, symbols, and units: report of the ASBMR histomorphometry nomenclature committee. J. Bone Miner. Res. 2, 595–610. doi: 10.1002/jbmr.5650020617

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazzaglia, U. E., Congiu, T., Marchese, M., and Dell'Orbo, C. (2010). The shape modulation of osteoblast-osteocyte transformation and its correlation with the fibrillar organization in secondary osteons. Cell Tissue Res. 340, 533–540. doi: 10.1007/s00441-010-0970-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazzaglia, U. E., Congiu, T., Sibilia, V., and Quacci, D. (2014). Osteoblast-osteocyte transformation. A SEM densitometric analysis of endosteal apposition in rabbit femur. J. Anat. 224, 132–141. doi: 10.1111/joa.12138

PubMed Abstract | CrossRef Full Text | Google Scholar

Polascik, T. J., and Mouraviev, V. (2008). Zoledronic acid in the management of metastatic bone disease. Therapeut. Clin. Risk Manage. 4:261. doi: 10.2147/TCRM.S2707

PubMed Abstract | CrossRef Full Text | Google Scholar

Poole, K. E., van Bezooijen, R. L., Loveridge, N., Hamersma, H., Papapoulos, S. E., Löwik, C. W., et al. (2005). Sclerostin is a delayed secreted product of osteocytes that inhibits bone formation. FASEB J. 19, 1842–1844. doi: 10.1096/fj.05-4221fje

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, L., Mak, A. T., Cheng, C., Hung, L., and Chan, K. (1999). Histomorphological study on pattern of fluid movement in cortical bone in goats. Anatom. Rec. 255, 380–387. doi: 10.1002/(SICI)1097-0185(19990801)255:4<380::AID-AR3>3.0.CO;2-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rejniak, K. A., and Anderson, A. R. A. (2011). Hybrid models of tumor growth. Wiley Interdiscipl. Rev. Syst. Biol. Med. 3, 115–125. doi: 10.1002/wsbm.102

PubMed Abstract | CrossRef Full Text | Google Scholar

Roudier, M. P., True, L. D., Higano, C. S., Vesselle, H., Ellis, W., Lange, P., et al. (2003). Phenotypic heterogeneity of end-stage prostate carcinoma metastatic to bone. Hum. Pathol. 34, 646–653. doi: 10.1016/S0046-8177(03)00190-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryser, M. D., Komarova, S. V., and Nigam, N. (2010). The cellular dynamics of bone remodeling: a mathematical model. SIAM J. Appl. Math. 70, 1899–1921. doi: 10.1137/090746094

CrossRef Full Text | Google Scholar

Ryser, M. D., Nigam, N., and Komarova, S. V. (2009). Mathematical modeling of spatio-temporal dynamics of a single bone multicellular unit. J. Bone Miner. Res. 24, 860–870. doi: 10.1359/jbmr.081229

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryser, M. D., Qu, Y., and Komarova, S. V. (2012). Osteoprotegerin in bone metastases: mathematical solution to the puzzle. PLoS Comput. Biol. 8:e1002703. doi: 10.1371/journal.pcbi.1002703

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapir-Koren, R., and Livshits, G. (2014). Osteocyte control of bone remodeling: is sclerostin a key molecular coordinator of the balanced bone resorption-formation cycles? Osteoporos. Int. 25, 2685–2700. doi: 10.1007/s00198-014-2808-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Skedros, J. G., Holmes, J. L., Vajda, E. G., and Bloebaum, R. D. (2005). Cement lines of secondary osteons in human bone are not mineral-deficient: new data in a historical perspective. Anatom. Rec. A Discov. Mol. Cell. Evol. Biol. 286A, 781–803. doi: 10.1002/ar.a.20214

PubMed Abstract | CrossRef Full Text | Google Scholar

Steck, R., and Knothe Tate, M. (2005). In silico stochastic network models that emulate the molecular sieving characteristics of bone. Ann. Biomed. Eng. 33, 87–94. doi: 10.1007/s10439-005-8966-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Stinson, J. C. (1975). The ailing mythical osteocyte. Med. Hypothes. 1, 186–190. doi: 10.1016/0306-9877(75)90049-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tauro, M., Shay, G., Sansil, S. S., Laghezza, A., Tortorella, P., Neuger, A. M., et al. (2017). Bone-seeking matrix metalloproteinase-2 inhibitors prevent bone metastatic breast cancer growth. Mol. Cancer Therapeut. 16, 494–505. doi: 10.1158/1535-7163.MCT-16-0315-T

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor-King, J. P., Basanta, D., Chapman, S. J., and Porter, M. A. (2017). Mean-field approach to evolving spatial networks, with an application to osteocyte network formation. Phys. Rev. E 96:012301. doi: 10.1103/PhysRevE.96.012301

PubMed Abstract | CrossRef Full Text | Google Scholar

van Oers, R. F. M., Klein-Nulend, J., and Bacabac, R. L. G. (2014). The osteocyte as an orchestrator of bone remodeling: an engineer's perspective. Clin. Rev. Bone Miner. Metab. 12, 2–13. doi: 10.1007/s12018-014-9154-9

CrossRef Full Text | Google Scholar

Wimpenny, D. I., and Moroz, A. (2007). On allosteric control model of bone turnover cycle containing osteocyte regulation loop. Biosystems 90, 295–308. doi: 10.1016/j.biosystems.2006.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, K., Xue, H., Dong, D., Zhu, S., Wang, J., Zhang, Q., et al. (2006). Identification of the proliferation/differentiation switch in the cellular network of multicellular organisms. PLoS Comput. Biol. 2:e145. doi: 10.1371/journal.pcbi.0020145

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, J., Piemontese, M., Thostenson, J. D., Weinstein, R. S., Manolagas, S. C., and O'Brien, C. A. (2014). Osteocyte-derived RANKL is a critical mediator of the increased bone resorption caused by dietary calcium deficiency. Bone 66, 146–154. doi: 10.1016/j.bone.2014.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Barragan-Adjemian, C., Ye, L., Kotha, S., Dallas, M., Lu, Y., et al. (2006). E11/gp38 selective expression in osteocytes: regulation by mechanical strain and role in dendrite elongation. Mol. Cell. Biol. 26, 4539–4552. doi: 10.1128/MCB.02120-05

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X. H.-F., Jin, X., Malladi, S., Zou, Y., Wen, Y. H., Brogi, E., et al. (2013). Selection of bone metastasis seeds by mesenchymal signals in the primary tumor stroma. Cell 154, 1060–1073. doi: 10.1016/j.cell.2013.07.036

PubMed Abstract | CrossRef Full Text | Google Scholar