A Mean-Field Approach to Evolving Spatial Networks, with an Application to Osteocyte Network Formation

TL;DR

This paper develops a mean-field framework for evolving spatial networks with node properties, deriving an integro-PDE for the local state degree distribution, and applies it to osteocyte network formation, revealing how differentiation rates affect network density and dendrite count.

Contribution

It introduces a novel mean-field approach for networks with evolving node properties and derives an integro-PDE for the local state degree distribution, validated through simulations.

Findings

Higher differentiation rates increase osteocyte density.

Increased differentiation leads to fewer dendrites per osteocyte.

The model aligns well with stochastic simulation results.

Abstract

We consider evolving networks in which each node can have various associated properties (a state) in addition to those that arise from network structure. For example, each node can have a spatial location and a velocity, or some more abstract internal property that describes something like social trait. Edges between nodes are created and destroyed, and new nodes enter the system. We introduce a "local state degree distribution" (LSDD) as the degree distribution at a particular point in state space. We then make a mean-field assumption and thereby derive an integro-partial differential equation that is satisfied by the LSDD. We perform numerical experiments and find good agreement between solutions of the integro-differential equation and the LSDD from stochastic simulations of the full model. To illustrate our theory, we apply it to a simple continuum model for osteocyte network formation within bones, with a view to understanding changes that may take place during cancer. Our results suggest that increased rates of differentiation lead to higher densities of osteocytes but with a lower number of dendrites. To help provide biological context, we also include an introduction to osteocytes, the formation of osteocyte networks, and the role of osteocytes in bona metastasis.

Equations

$$u_k(t,\mathbf{s})=p_k(t|\mathbf{s})f(t,\mathbf{s}),$$

$$f(t,\mathbf{s})=\sum _{k=0}^{\infty }{u}_k(t,\mathbf{s}),p_k(t|\mathbf{s})=\frac{u_k(t,\mathbf{s})}{{\sum }_{k=0}^{\infty }{u}_k(t,\mathbf{s})}.$$

$$P_k(t)=\frac{{\int }_{\text{mathds}S}{u}_k(t,\mathbf{s})\text{d}\mathbf{s}}{{\sum }_{k=0}^{\infty }{\int }_{\text{mathds}S}{u}_k(t,\mathbf{s})\text{d}\mathbf{s}}.$$

$${\dot{\mathbf{x}}}_i={\mathbf{v}}_i,m{\dot{\mathbf{v}}}_i=-\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{\nabla }_{{\mathbf{x}}_i}\Phi ({\mathbf{x}}_i-{\mathbf{x}}_j).$$

$$\text{d}{\mathbf{X}}_i=\mathbf{\mu }({\mathbf{X}}_i)\text{d}t+\sigma \text{d}{\mathbf{W}}_t,$$

$$\sum _{n=0}^{\infty }\left\{\sum _{K_{1,n}}\left[{\int }_{\text{mathds}{S}^n}{\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)\text{d}{\vec{s}}_n\right]\right\}=1,$$

$$\sum _{K_{a,b}}\equiv \sum _{k_a=0}^{\infty }\sum _{k_{a+1}=0}^{\infty }\cdots \sum _{k_b=0}^{\infty }.$$

$$\left(\frac{\partial }{\partial t}+{\mathcal{L}}^{(n)}\right){\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)=\sum _{i=1}^n\sum _{j=i+1}^n\left(\mathcal{C}({\mathbf{s}}_i,k_i-1,{\mathbf{s}}_j,k_j-1){\mathcal{F}}_n^{{\vec{k}}_{n,-}^{i,j}}(t,{\vec{s}}_n)-\mathcal{C}({\mathbf{s}}_i,k_i,{\mathbf{s}}_j,k_j){\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)\right)$$

$$+\sum _{i=1}^n\frac{1}{n}{\delta }_{k_i,0}\mathcal{J}\mathcal{P}({\mathbf{s}}_i){\mathcal{F}}_{n-1}^{{\vec{k}}_n^{i-}}(t,{\vec{s}}_n^{i-})-\mathcal{J}{\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n).$$

$${\mathcal{L}}^{(n)}{\mathcal{F}}_n^{{\vec{k}}_n}=\sum _{i=1}^n\left({\mathbf{v}}_i\cdot {\nabla }_{{\mathbf{x}}_i}+\frac{\mathbf{F}({\vec{x}}_n)}{m}\cdot {\nabla }_{{\mathbf{v}}_i}\right){\mathcal{F}}_n^{{\vec{k}}_n},$$

$$\mathbf{F}({\vec{x}}_n)=-\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{\nabla }_{{\mathbf{x}}_i}\Phi ({\mathbf{x}}_i-{\mathbf{x}}_j).$$

$${\mathcal{L}}^{(n)}{\mathcal{F}}_n^{{\vec{k}}_n}=\sum _{i=1}^n{\nabla }_{{\mathbf{x}}_i}\cdot \left(\mathbf{\mu }({\mathbf{x}}_i)-\frac{{\sigma }^2}{2}{\nabla }_{{\mathbf{x}}_i}\right){\mathcal{F}}_n^{{\vec{k}}_n}.$$

$$u_{k_1}(t,s_1)=\sum _{n=0}^{\infty }n\sum _{K_{2,n}}{\int }_{\text{mathds}{S}^{n-1}}{\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)\text{d}{\vec{s}}_n^{(2)}=\sum _{n=0}^{\infty }n\sum _{k_2=0}^{\infty }\cdots \sum _{k_n=0}^{\infty }{\int }_{\text{mathds}{S}^{n-1}}{\mathcal{F}}_n^{{\vec{k}}_n}(t,{\mathbf{s}}_1,\dots ,{\mathbf{s}}_n)\text{d}{\mathbf{s}}_2\dots \text{d}{\mathbf{s}}_n,$$

$$\sum _{n=0}^{\infty }n\sum _{K_{2,n}}{\int }_{\text{mathds}{S}^{n-1}}\frac{\partial {\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)}{\partial t}\text{d}{\vec{s}}_n^{(2)}=\frac{\partial u_{k_1}(t,{\mathbf{s}}_1)}{\partial t}.$$

$$\sum _{n=0}^{\infty }n\sum _{K_{2,n}}{\int }_{\text{mathds}{S}^{n-1}}{\mathcal{L}}^{(n)}{\mathcal{F}}_n^{{\vec{k}}_n}\text{d}{\vec{s}}_n^{(2)}.$$

$${\mathcal{L}}^{(1)}{u}_{k_1}={\nabla }_{{\mathbf{x}}_1}\cdot \left(\mathbf{\mu }({\mathbf{x}}_1)-\frac{{\sigma }^2}{2}{\nabla }_{{\mathbf{x}}_1}\right)u_{k_1}.$$

$$u_{k_1,k_2}^{(2)}(t,{\mathbf{s}}_1,{\mathbf{s}}_2)=\sum _{n=0}^{\infty }n(n-1)\sum _{K_{3,n}}{\int }_{\text{mathds}{S}^{n-2}}{\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)\text{d}{\vec{s}}_n^{(3)}.$$

$$u_{k_1,k_2}^{(2)}(t,{\mathbf{s}}_1,{\mathbf{s}}_2)\approx u_{k_1}(t,{\mathbf{s}}_1)u_{k_2}(t,{\mathbf{s}}_2).$$

$${\bar{\mathcal{L}}}^{(1)}{u}_{k_1}={\mathbf{v}}_1\cdot {\nabla }_{{\mathbf{x}}_1}{u}_{k_1}-\mathcal{B}(f,u_{k_1}),$$

$$\mathcal{B}(f,u_{k_1})=\left({\nabla }_{{\mathbf{x}}_1}\Phi \ast {\int }_{\text{mathds}{R}^d}f\text{d}{\mathbf{v}}_1\right)\cdot \left({\nabla }_{{\mathbf{v}}_1}{u}_{k_1}\right),$$

$$={\int }_{\text{mathds}S}\sum _{k_2=0}^{\infty }\left(\mathcal{C}({\mathbf{s}}_1,k_1-1,{\mathbf{s}}_2,k_2-1)u_{k_1-1,k_2-1}^{(2)}(t,{\mathbf{s}}_1,{\mathbf{s}}_2)-\mathcal{C}({\mathbf{s}}_1,k_1,{\mathbf{s}}_2,k_2)u_{k_1,k_2}^{(2)}(t,{\mathbf{s}}_1,{\mathbf{s}}_2)\right)\text{d}{\mathbf{s}}_2+\mathcal{J}\mathcal{P}({\mathbf{s}}_1){\delta }_{k_1,0}.$$

$$\left[\frac{\partial }{\partial t}+{\mathcal{L}}^{(1)}\right]u_{k_1}(t,{\mathbf{s}}_1)$$

$$-\left({\int }_{\text{mathds}S}\sum _{k_2=0}^{\infty }\mathcal{C}({\mathbf{s}}_1,k_1,{\mathbf{s}}_2,k_2)u_{k_2}(t,{\mathbf{s}}_2)\text{d}{\mathbf{s}}_2\right)u_{k_1}(t,{\mathbf{s}}_1)+\mathcal{J}\mathcal{P}({\mathbf{s}}_1){\delta }_{k_1,0},$$

$$\sum _{n=1}^{\infty }\left\{\sum _{A_n\in S_n}\left[{\int }_{\text{mathds}{S}^n}{\mathcal{F}}_n^{A_n}(t,{\vec{s}}_n)\text{d}{\vec{s}}_n\right]\right\}=1,$$

$$S_{n}=\left\{A_{n}:\begin{array}[]{ll}(A_{n})_{ij}\in\mathds{N}_{0}\,,&1\leq i<j\leq n\\[5.69054pt] (A_{n})_{ij}=0\,,&\mbox{ otherwise}\end{array}\right\}\,,$$

$$\sum _{A_n\in S_n}\equiv \sum _{(A_n{)}_{12}=0}^{\infty }\cdots \sum _{(A_n{)}_{1n}=0}^{\infty }\sum _{(A_n{)}_{23}=0}^{\infty }\cdots \sum _{(A_n{)}_{n-1,n}=0}^{\infty }.$$

$$k_i=\sum _{j=i+1}^n(A_n{)}_{ij}+\sum _{j=1}^{i-1}(A_n{)}_{ji}.$$

$${\mathcal{F}}_n^{{\vec{k}}_n}(t,{\vec{s}}_n)$$

$$\left(\frac{\partial }{\partial t}+{\mathcal{L}}^{(n)}\right){\mathcal{F}}_n^{A_n}(t,{\vec{s}}_n)=\sum _{i=1}^n\sum _{j=i+1}^n\mathcal{C}({\mathbf{s}}_i,k_i-1,{\mathbf{s}}_j,k_j-1){\mathcal{F}}_n^{A_{n,-}^{ij}}(t,{\vec{s}}_n)-\mathcal{C}({\mathbf{s}}_i,k_i,{\mathbf{s}}_j,k_j){\mathcal{F}}_n^{A_n}(t,{\vec{s}}_n)$$

$$+\sum _{i=1}^n\sum _{j=i+1}^n[(A_n{)}_{ij}+1]\mathcal{D}({\mathbf{s}}_i,k_i+1,{\mathbf{s}}_j,k_j+1){\mathcal{F}}_n^{A_{n,+}^{ij}}(t,{\vec{s}}_n)-(A_n{)}_{ij}\mathcal{D}({\mathbf{s}}_i,k_i,{\mathbf{s}}_j,k_j){\mathcal{F}}_n^{A_n}(t,{\vec{s}}_n)$$

$$+\frac{1}{n}\sum _{i=1}^n\left[\prod _{j=1}^{i-1}\delta (0,(A_n{)}_{ji})\right]\left[\prod _{j=i+1}^n\delta (0,(A_n{)}_{ij})\right]\mathcal{J}\mathcal{P}({\mathbf{s}}_i){\mathcal{F}}_{n-1}^{A_n^{i-}}(t,{\vec{s}}_n^{i-})-\mathcal{J}{\mathcal{F}}_n^{A_n}(t,{\vec{s}}_n),$$

$$(A_{n,\pm}^{ij})_{lm}=\left\{\begin{array}[]{ll}(A_{n})_{lm}\pm 1\,,&\text{if $(i,j)=(l,m)$}\,,\\ (A_{n})_{lm}\,,&\text{otherwise\,.}\end{array}\right.$$

Full text

I Introduction

Networks, in which entities (nodes'') interact with each other via edges'', are a useful representation of complex systems Newman_2010 . They have often been very helpful for formulating and answering questions in biology, sociology, engineering, and numerous other fields. Because the present work is motivated by a biological application, let's consider a few examples from biology. Many systems—such as blood vasculature Jain_1988 , leaf venation Blonder_2011 , and fungi heaton2010 ; shl2015 —can be treated as biological transportation networks, in which edges carry resources and nodes operate as junctions. Some of these studies exploit ideas from fluid mechanics and energy minimization to investigate flow through various media Liesbet_2013 ; Hu_2013 ; Haskovec_2015 ; heaton2010 . To give another type of example, in evolutionary game theory, a node can represent a biological agent, and edges indicate interactions in a ``game'' between those agents Szabo_2007 . Applications range from behavioral ecology Smith_1982 to investigating tumor heterogeneity in cancer Archetti_2013 . Changes in node fitness can depend on, for example, a node's phenotype and its immediate neighbors Ohtsuki_2006 .

Many of the above examples involve ``spatial networks'' Barthelemy_2011 , and spatial constraints can exert significant influence (directly and/or indirectly) on both network structure and function. In the aforementioned examples, the networks are embedded in space, and one thus can assign physical locations to the nodes and edges. This is clearly important when considering dynamical processes on those networks Durrett_1994 .

Some spatial networks grow in time as they form: new nodes and edges can join a spatial network, and the spatial domain can expand. For example, cities often grow outwards or arise when borders from multiple settlements coalesce Gallotti_2014 ; Louf_2014 ; barth2017 . Fungi, which are living networks, expand to reach nutrients, and such growth induces flows of mass heaton2010 ; shl2015 .

In the present paper, we propose a framework for describing evolving spatial networks. We then apply this framework to examine the formation of osteocyte networks in bone (see Fig. 8 for a schematic). During bone formation, cells called osteoblasts secrete bone matrix and differentiate into cells called osteocytes, and the ensuing growth process results in a network of connected cells that communicate by chemical diffusion via gap junctions Buenzli_2015 ; Dallas_2010 ; Franz-Odendaal_2006 ; Ishihara_2012 . In Section VII, we develop and analyze a model for this process, with motivation of using network analysis to study bone cancer. In pathological bone, the highly regulated bone-remodeling signaling pathway is disrupted, and it may be possible to gain insight into the nature of this disruption using tools from network analysis.

There are myriad models of network formation Newman_2010 . There are at least three possible ways of formulating such a model: (i) all nodes and edges are created simultaneously with a single algorithmic step (e.g., the standard Erdős–Rényi (ER) random graph $G(n,p)$ Erdos_1960 ; Newman_2010 and standard random geometric graphs penrosergg ); (ii) nodes and edges have an implicit order of creation but time is not considered explicitly (e.g., in some preferential-attachment models Price_1986 ; Moore_2006 ); or (iii) nodes and edges have an order of creation and time is considered explicitly (e.g., in some preferential-attachment models Krapivsky_2001 ; Krapivsky_2001_b and in adaptive network models thilo-adapt2 ). Because we want to incorporate time explicitly, we will consider spatial networks in category (iii) in the present paper. See also the recent work by Zuev et al. on geometric preferential-attachment models bianconi2015 .

When studying a model in category (iii), it is common to employ kinetic approximations Krapivsky_2010 ; Krapivsky_2001 ; Krapivsky_2001_b ; Marceau_2010 . Such approximations often allow one to construct a ``master equation'' to obtain approximate and/or asymptotic expressions for quantities such as degree distributions, component sizes, and cycle sizes Krapivsky_2010 . These equations can take the form of an ordinary differential equation (ODE), partial differential equation (PDE), or other continuous model. By carefully constructing a general state space of the system, we derive an extension to the ODE master equations given in Krapivsky_2010 ; Krapivsky_2001 ; Krapivsky_2001_b ; Marceau_2010 to obtain a master equation in the form of an integro-partial differential equation (IPDE) that incorporates this state space.

Our work illustrates how to use a master-equation approach to study spatial networks. As an illustration of its potential, we examine degree distributions in a very general model of evolving spatial networks and use our results to gain insights on osteocyte network formation in bone. An important benefit of using an explicit time-dependent kinetic approach is that it allows one to incorporate nodes that move through a state space (e.g., particles that diffuse). Note that we will often use the terms nodes'' and particles'' interchangeably (depending on the context). When considering spatially-embedded network models, the state space corresponds to each node being located in a copy of the physical space. From a mathematical viewpoint, our approach is reminiscent of some models of social networks Boguna_2003 ; Wong_2006 , for which the state space corresponds to a latent social space described by some internal parameters; however, these models have no time dependence. See also the recent paper barre2017 on macroscopic descriptions for particle interactions mediated by time-dependent networks.

The remainder of our paper is organized as follows. In Section II, we define the concept of a ``local state degree distribution'' (LSDD), which encapsulates the degree distribution of a network local to a point in state space. In Section III, we give a description of our model of evolving spatial networks. In Section IV, we derive an equation for the LSDD when edges are created but cannot be destroyed. In Section V, we incorporate edge destruction into our model. Our derivations require the use of approximations, so we use numerical simulations to explore agreement and discrepancies between theory and our derived equations in Section VI. In Section VII, we use our model for evolving spatial networks to formulate a model for osteocyte network formation. We conclude and discuss future directions in Section VIII. In appendices, we give additional details on derivations and numerical implementations.

II A Local State Degree Distribution

We consider networks in which each node (i.e., particle) has some number of associated properties in addition to those, such as degree distribution, that arise from network structure. For example, each node may have a spatial location and a velocity, or it may have some more abstract internal (``latent'') property describing, for example, some social trait. We collect these properties together into a state vector $\boldsymbol{s}$, which belongs to a state space $\mathds{S}$.

The density $f(t,\boldsymbol{s})$ of particles in the state space gives the expected number of particles with state $\boldsymbol{s}$ 111More precisely, $f(t,\boldsymbol{s})\text{d}\boldsymbol{s}$ gives the expected number of particles that have states lying in the volume element $\text{d}\boldsymbol{s}$ centered on $\boldsymbol{s}$.. However, a common way to study the properties of models of network formation is to examine the degree distribution Newman_2010 . In this paper, we combine these ideas to consider what we call a local state degree distribution (LSDD) $u_{k}(t,\boldsymbol{s})$, which gives the expected number of particles of degree $k$ at time $t$ with state vector $\boldsymbol{s}$. One can write the LSDD as

[TABLE]

where $p_{k}(t\,|\,\boldsymbol{s})$ is the conditional probability that a node at time $t$ has degree $k$, given that its state is $\boldsymbol{s}$. We can recover both $f$ and $p_{k}$ from $u_{k}$ via

[TABLE]

The degree distribution of the whole network is given by

[TABLE]

III Model of Evolving Spatial Networks

We now present a model for evolving spatial networks. Our model has three tunable features. First, edges can be created and deleted. Second, new nodes can be created (but we do not allow node deletion). Finally, we specify a model (possibly depending on network structure) for the evolution of the state $\boldsymbol{s}$ of each node.

Suppose at time $t$ that there are $N(t)$ nodes with state vectors $\boldsymbol{s}_{i}$ and degrees $k_{i}$ (with $i\in\{1,\ldots,N\}$). We also suppose that new edges are created between each pair of nodes as independent Poisson processes, where $\mathcal{C}(\boldsymbol{s}_{i},k_{i},\boldsymbol{s}_{j},k_{j})$ is the rate of edge creation between node $i$ and node $j$, so that the probability of an edge being created between node $i$ and node $j$ in time $t$ to $t+\text{d}t$ is $\mathcal{C}(\boldsymbol{s}_{i},k_{i},\boldsymbol{s}_{j},k_{j})\text{d}t$. We also suppose that $\mathcal{C}$ depends on the states and degrees of the two nodes $i$ and $j$, but that it does not depend on other properties of the network (such as, for example, whether an edge between node $i$ and node $j$ already exists). Thus, our model allows multiedges (i.e., multiple edges between two distinct nodes).

Similarly, we suppose that existing nodes are deleted as independent Poisson processes, where $\mathcal{D}(\boldsymbol{s}_{i},k_{i},\boldsymbol{s}_{j},k_{j})$ is the rate of edge deletion (per edge) between node $i$ and node $j$. We suppose also that new particles, which have degree [math], arrive randomly as a Poisson process with constant rate $\mathcal{J}$, and we assign to them a state drawn at random from the probability distribution $\mathcal{P}$.

The final component of our model is the equation of motion of the particles in the state space. Many possible models are available (both deterministic and stochastic), and we want to keep our presentation as general as possible. Nevertheless, it is useful to have a model in mind to fix ideas, and we thus present several examples.

The simplest case is for each node to have a constant state vector. Our first nontrivial example consists of identical particles of mass $m>0$ that follow Newton's equations of motion with a smooth pairwise potential $\Phi$. In that case, $\boldsymbol{s}_{i}(t)=(\boldsymbol{x}_{i}(t),\boldsymbol{v}_{i}(t))\in\mathds{S}\subseteq\mathds{R}^{2d}$ and

[TABLE]

Variations of these equations have been used to consider collective motion such as swarming Degond_2004 ; Gallagher_2013 ; Illner_1987 . Our second example consists of the stochastic position jump process in which $\boldsymbol{s}_{i}(t)=\boldsymbol{X}_{i}(t)\in\mathds{S}\subset\mathds{R}^{d}$ and 222Following standard practice, we use capital letters for random variables and lower-case letters for realizations of these variables.

[TABLE]

where $\boldsymbol{W}_{t}$ is standard Brownian motion (i.e., a Wiener process) Oksendal_2014 . Other stochastic examples include velocity-jump processes Othmer_1988 and fractional diffusion processes Metzler_2000 .

In each of the above examples, the motion in state space is independent of network structure. Of course, it is also possible to imagine scenarios in which the motion depends on node degree or other structural features.

We summarize the model events and state update in Table 1, and we illustrate them in Fig. 1.

Our model treats all nodes identically, although our methodology can be extended to heterogeneous classes of nodes or heterogeneous classes of edges. (Both of these generalizations are examples of multilayer networks Kivela_2014 .) To consider such cases, one can modify kinetic methodology from Bisi_2006 ; Burgers_1961 .

IV Model Analysis: No Edge Deletion

We begin by considering the case in which there is no edge deletion. That is $\mathcal{D}(\boldsymbol{s}_{i},k_{i},\boldsymbol{s}_{j},k_{j})=0$ for all $\boldsymbol{s}_{i},k_{i},\boldsymbol{s}_{j},k_{j}$. In the next subsection, we give a set of hierarchical master Fokker–Planck (FP) equations for the probability distribution of the state of the system. Because it is not pragmatic to work in this high-dimensional space, in Section IV.2, we reduce the dimension using mean-field arguments from kinetic theory.

IV.1 Fokker–Planck Equation

We define $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$ to be the probability that a network has $n$ nodes with degree sequence $\vec{k}_{n}=\{k_{1},\dots,k_{n}\}$ and state vectors $\vec{s}_{n}=\{\boldsymbol{s}_{1},\dots,\boldsymbol{s}_{n}\}$. Note the normalization

[TABLE]

where $K_{a,b}$ is shorthand for summing over all possible degrees for nodes $i\in\{a,a+1,\dots,b\}$. That is,

[TABLE]

We are not considering edge deletion and we have assumed that edge creation and state-space motion depend only on node state and degree, so it is possible to write down a closed equation for $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$. In contrast, when we do consider edge deletion in Section V, it will not be enough simply to keep track of node degrees. We will need the full adjacency matrix.

Because the probability density function $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$ depends on the number of particles in the system (which changes when new nodes are created), we obtain a hierarchy of Fokker–Planck equations as in Chen_2014 . By considering a small time step from $t$ to $t+\text{d}t$ and partitioning over the events that can occur, we obtain

[TABLE]

The operator on the left-hand side describes the evolution of particles in state space, and it depends on the particular model that one chooses for this evolution. For example, if nodes move according to equation (3), one has the Liouville flux term

[TABLE]

where $\vec{x}_{n}=(\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})$ and

[TABLE]

Alternatively, if $\boldsymbol{s}_{i}(t)$ evolves according to the stochastic differential equation (SDE) (4), one has the Kolmogorov forward operator

[TABLE]

If the states of nodes are time-independent, then $\mathcal{L}^{(n)}\equiv 0$.

The first term in parentheses on the right-hand side of equation (7) corresponds to edge-creation events between each pair, $i$ and $j$, of nodes. The positive term corresponds to gaining a network with degree sequence $\vec{k}_{n}=\{k_{1},\dots,k_{n}\}$ from a network with degree sequence $\vec{k}_{n,-}^{i,j}=\{k_{1},\dots,k_{i}-1,\ldots,k_{j}-1,\dots,k_{n}\}$ by adding an edge between nodes $i$ and $j$. The negative term corresponds to losing a network of degree sequence $\vec{k}_{n}$ (as it changes to a network of degree sequence $\vec{k}_{n,+}^{i,j}=\{k_{1},\dots,k_{i}+1,\dots,k_{j}+1,\dots,k_{n}\}$) by adding an edge between nodes $i$ and $j$.

The second term on the right-hand side of equation (7) corresponds gaining a network with degree sequence $\{k_{1},\dots,k_{i-1},0,k_{i+1},\dots,k_{n}\}$ from a network with degree sequence $\vec{k}_{n}^{i-}=\{k_{1},\dots,k_{i-1},k_{i+1},\dots,k_{n-1}\}$ by adding a new node (of degree [math]) to the system. The Kronecker delta $\delta_{k_{i},0}$ ensures that this term is present only when $k_{i}=0$, corresponding to the new node having degree [math]. One draws the new state $\boldsymbol{s}_{i}$ from the probability distribution with density function $\mathcal{P}$. (The state vector of the existing nodes is $\vec{s}_{n}^{i-}=(\boldsymbol{s}_{1},\dots,\boldsymbol{s}_{i-1},\boldsymbol{s}_{i+1},\dots,\boldsymbol{s}_{n})$.) We assign the label of the new node uniformly at random from the set $\{1,\ldots,n\}$ (rather than assigning it to be the last node $n$) to ensure that particles are indistinguishable.

The final term on the right-hand side of equation (7) describes the loss of a network with degree sequence $\vec{k}_{n}$ because of the addition of a new node. (One thereby obtains a network with degree sequence $\{k_{1},\dots,k_{i-1},0,k_{i},\dots,k_{n}\}$ for some $i$.)

To save us from writing down separate equations for each case in which $k_{i}=0$ for some $i$ (because it is impossible to arrive at a state in which a node has degree [math] by adding an edge to a state in which it has degree $-1$), we use the convention that $\mathcal{F}_{n}^{\vec{k}_{n}}(t,\vec{s}_{n})=0$ if $n<0$ or $k_{i}<0$ for all $i\in\{0,1,\dots,n\}$. We also suppose that particles are indistinguishable initially, so the initial condition $\mathcal{F}_{n}^{\,\vec{k}_{n}}(0,\vec{s}_{n})$ is invariant to index permutation. Equation (7) then ensures that this is true for all $t$.

It is not feasible to solve equation (7) analytically (except perhaps when the node state vectors are uncorrelated), and it is not practical to solve it numerically due to the high dimension of the domain. In Section IV.2, we reduce the dimension of the equation (7) using mean-field approaches from kinetic theory.

IV.2 Low-Dimensional Approximation

To derive our low-dimensional approximation, we adapt methods from kinetic theory Carillo_2009 ; Degond_2004 . We keep the presentation brief in this subsection; we give more details in Appendix B.

A common approach in kinetic theory is to average over the states of particles $2$ to $n$ to find an equation for the marginal distribution function of the first particle Cercignani_1988 ; Harris_1971 (the so-called ``1-particle distribution function''). Because particles are indistinguishable, multiplying by $n$ gives the probability of finding any particle in a given state. Here we adopt the same approach, and we average over the states and degrees of particles $2$ to $n$. Because the number of particles itself can vary, we also need to average over this quantity. The resulting 1-particle distribution function is exactly the previously-defined LSDD $u_{k}(t,s)$. Specifically,

[TABLE]

where we introduce the shorthand notation $\text{d}\vec{s}_{n}^{\,(\mu)}=\text{d}\boldsymbol{s}_{\mu}\dots\text{d}\boldsymbol{s}_{n}$ for $\mu\in\{1,\dots,n\}$.

To find the equation satisfied by $u_{k}(t,s_{1})$, we apply the same summation and integration to the Fokker–Planck equation (7).

Because the summation and integration commutes with the time derivative, it follows for the first term on the left-hand side (LHS) of equation (7) that

[TABLE]

For the next term on the LHS, we need to evaluate

[TABLE]

If there are no interactions between particles in the state space, this term evaluates to $\mathcal{L}^{(1)}u_{k_{1}}$. Thus, for example, if $\mathcal{L}^{(n)}$ is given by (9), then this term is

[TABLE]

When there are pairwise interactions between particles in the state space, for each interacting pair, one can perform the integration over all other particles. Consequently, after relabelling, one can write (12) in terms of the 2-particle LSDD

[TABLE]

We do not have a closed equation for $u_{k_{1}}$, but the first in a series of equations (the BBGKY hierarchy) for the 1-particle, 2-particle, 3-particle, etc. LSDDs. In this case, we make the common mean-field closure assumption that

[TABLE]

For example, if $\mathcal{L}^{(n)}$ is given by equation (8), then equation (12) becomes

[TABLE]

where $\mathcal{B}(f,u_{k_{1}})$ is the mean-field approximation

[TABLE]

where $*$ represents the convolution operator and the function $f$ is given in equation (1).

We now apply the same integration and summation to the right-hand side (RHS) of equation (7). We give a detailed derivation in Appendix B. Here we simply note that by relabelling particles and again using the fact that $\mathcal{F}_{n}^{\vec{k}_{n}}(t,\vec{s}_{n})$ is invariant with respect to index permutation, we find that

[TABLE]

We again need to use the mean-field closure assumption (15) to write the 2-particle LSDD in terms of the 1-particle LSDD. This gives the final closed mean-field equation for the 1-particle LSDD in the absence of edge deletion:

[TABLE]

where $u_{k_{1}}\equiv 0$ if $k_{1}<0$ by convention. When $\mathcal{C}$ is a constant and $\mathds{S}$ is a point, equation (IV.2) reduces to the master equations given in Ref. Krapivsky_2010 . The quadratic terms in equation (IV.2) are analogous to the mean-field term in the Vlaslov equation, where a test particle feels the effect of a ``cloud'' of points Cercignani_1988 ; Harris_1971 ; Gallagher_2013 .

V Derivation of the Model: Edge Deletion

The state space $\{(\vec{s}_{n},\vec{k}_{n}):n>0\}$ that we used in Section IV is not sufficient when we allow edge deletion. With edge deletion, it is crucial to know whether an edge exists between each pair of nodes, so we must consider the underlying adjacency matrix. For undirected networks with multiedges, the adjacency matrix $A_{n}=A_{n}^{T}$ has entries $(A_{n})_{i,j}$ for $i\neq j$ and $i,j\in\{1,\dots,n\}$, where $(A_{n})_{i,j}\in\mathds{N}_{0}$ gives the number of edges between nodes $i$ and $j$.

Because we consider probability distributions over $A_{n}$, the most efficient representation is to restrict attention to the independent entries of $A_{n}$. We thus change convention slightly and set $(A_{n})_{i,j}=0$ for $i\geq j$, and we will retain the term ``adjacency matrix'' to indicate the resulting matrix.

Let $\mathscr{F}_{n}^{A_{n}}(t,\vec{s}_{n})$ denote the probability that a network has adjacency matrix $A_{n}$ and $n$ nodes with state vectors $\vec{s}_{n}=\{\boldsymbol{s}_{1},\dots,\boldsymbol{s}_{n}\}$. The normalization condition is

[TABLE]

where

[TABLE]

so that

[TABLE]

For fixed $n$, one can calculate the degree $k_{i}$ of node $i$ from the adjacency matrix using

[TABLE]

Therefore, we can relate the distributions $\mathcal{F}$ and $\mathscr{F}$ via

[TABLE]

where $\delta(a,b)$ is the Kronecker delta (which is usually written as $\delta_{ab}$). When we include edge deletion, the hierarchical Fokker–Planck equation (7) becomes

[TABLE]

where

[TABLE]

The first term on the RHS of equation (23) corresponds to edge-creation events between nodes $i$ and $j$ as before (see Section IV). The first part of the second term corresponds to gaining a network with adjacency matrix $A_{n}$ from a network with adjacency matrix $A_{n,+}^{ij}$ by deleting an edge between $i$ and $j$. Note that $\mathcal{D}$ is the rate of deletion per edge, so we multiply by the number of edges (which is equal to $(A_{n})_{ij}+1$) between $i$ and $j$. The second part of this term corresponds to losing a network with adjacency matrix $A_{n}$ by deleting an edge between $i$ and $j$ (to produce a network with adjacency matrix $A_{n,-}^{ij}$).

The third term on the RHS of equation (23) corresponds to gaining a network with adjacency matrix

[TABLE]

from an adjacency matrix

[TABLE]

by relabeling nodes $j\rightarrow j+1$ for $j\geq i$ and adding a new unconnected node with label $i$ (which we choose uniformly at random from the set $\{1,\ldots,n\}$). The Kronecker $\delta$ ensures that this term is present only when $(A_{n})_{ij}=(A_{n})_{ji}=0$. As in Section IV, the uniformly random choice of the label for the new node ensures that $\mathscr{F}_{n}^{A_{n}}(t,\vec{s}_{n})$ is invariant with respect to index permutations.

We give the derivation of a reduced equation for $u_{k_{1}}(t,\boldsymbol{s}_{1})$ in Appendix C. When we include edge deletion, we find that we need an additional closure assumption in addition to the mean-field approximation (15); see Appendix C for details. Our low-dimensional approximation to (23) is

[TABLE]

VI Numerical Examples

We now carry out Monte-Carlo simulations of our stochastic network evolution process to illustrate the validity of equations (IV.2) and (25) for sufficiently large networks.

First, in Section VI.1, we investigate numerically the validity of the mean-field assumption for our model when the state space is a point. Second, in Section VI.2, we use an example scenario to demonstrate that a numerical solution of equations (IV.2) and (25) matches well with a full Monte-Carlo simulation of the underlying process. Third, in Section VI.3, we adapt the example scenario to consider the convergence of a network's degree distribution in the limit of large networks. Fourth, Section VI.4, we show that one can use our kinetic approximation as an alternative to some one-step network creation models proposed by Boguñá et al. Boguna_2003 . Finally, in Section VI.5, we briefly consider some further approximations that one can make to equations (IV.2) and (25).

VI.1 Accuracy of Our Mean-Field Assumption

In the absence of a network structure, the validity of the mean-field approximation depends on the choice of $\mathcal{L}^{(n)}$, and it has been studied widely Cercignani_1988 ; Harris_1971 ; Gallagher_2013 ; Bruna_2012 ; Degond_2004 . In the present paper, we focus on evolving network structure. In our exploration of the validity of the mean-field approximation in the network evolution model, we suppose that there is no state dependence in either the node creation rate $\mathcal{C}$ or the node deletion rate $\mathcal{D}$.

The 2-particle degree distribution $P^{(2)}_{k_{1},k_{2}}$ is the probability that two nodes selected uniformly at random without replacement have degrees $k_{1}$ and $k_{2}$. In our mean-field closure, we approximate this quantity by the product $P_{k_{1}}P_{k_{2}}$, where $P_{k_{1}}$ is the probability that a single node selected uniformly at random has degree $k_{1}$ (and $P_{k_{2}}$ is defined analogously). In Fig. 2, we compare the empirical distributions $P_{k_{1},k_{2}}^{(2)}(t)$ and $P_{k_{1}}(t)P_{k_{2}}(t)$ generated from 100 realizations of Algorithm 1 using 125 particles (and no node creation).

We see that our mean-field approximation does well on this example, and the main error occurs when $k_{1}=k_{2}$, for which the product $P_{k_{1}}(t)P_{k_{2}}(t)$ of 1-particle distributions is slightly larger than the 2-particle distribution $P_{k_{1},k_{2}}^{(2)}(t)$. This discrepancy arises because one cannot select the same node twice when evaluating the correlation function $P_{k_{1},k_{2}}^{(2)}(t)$, so the probability of finding two nodes with the same degree is lower than that estimated by $P_{k_{1}}(t)P_{k_{2}}(t)$ (which corresponds to choosing two nodes uniformly randomly with replacement). The difference should therefore tend to [math] as $1/n$ as the number $n$ of nodes becomes infinite.

One way to evaluate the difference between two probability distributions is the Kolmogorov–Smirnov (KS) test Peacock_1983 , which gives the probability $\rho_{\text{KS}}$ that one rejects the hypothesis that the two distributions are equal. For $P_{k_{1},k_{2}}^{(2)}(t)$ and $P_{k_{1}}(t)P_{k_{2}}(t)$, we find that $\rho_{\text{KS}}\approx 2.0\times 10^{-3}$.

VI.2 Example Scenario: Local State Degree Distribution

The example in Section VI.1 was particularly simple, as it focused only on the network aspect of the model. We now want to compare Monte-Carlo simulations of a much more complicated evolving spatial network with a numerical solution of the reduced equation (25). We select our new example to illustrate and evaluate all of the model components described in Table 1. It does not represent any particular physical or biological process.

Let's consider noninteracting point particles contained in the unit square, so one can describe the state of each particle by its position vector $\boldsymbol{s}_{i}=(x_{i},y_{i})\in[0,1]^{2}$. We suppose that these positions evolve according to the SDEs Oksendal_2014

[TABLE]

where as before we use capital letters $\boldsymbol{S}_{i}=(X_{i},Y_{i})$ to distinguish random variables from the values that they take. We assume that the drift coefficient $\mu>0$ and volatility coefficient $\sigma>0$ are constant (corresponding to a constant diffusion coefficient $\sigma^{2}/2$). To initialise, we place $1000$ particles with degree [math] uniformly at random in the rectangle $(X_{i},Y_{i})\in[0,1/10]\times[0,1]$. We impose reflective boundary conditions at $x=0$ and $x=1$, and we impose periodic boundary conditions at $y=0$ and $y=1$. To generate some spatial heterogeneity, we suppose that the rate of edge creation between nodes depends both on the distance between nodes and on the spatial coordinates of each node. We take

[TABLE]

In contrast, we suppose that the rate of edge deletion per edge between nodes $i$ and $j$ depends on the degrees of nodes $i$ and $j$ but is independent of position. Specifically, we take

[TABLE]

We introduce new nodes of degree [math] at a rate $\mathcal{J}$ uniformly at random in the rectangle $(x,y)\in[0,1/10]\times[0,1]$. Therefore,

[TABLE]

We give a schematic illustration of these processes in Fig. 3. We simulate the system using Algorithm 1 until final time $T_{\text{end}}=1/2$.

Having defined the stochastic process that we are simulating, we now turn to the reduced model (25). Because $\mathcal{C}(\boldsymbol{s}_{1},k_{1},\boldsymbol{s}_{2},k_{2})$ and $\mathcal{D}(\boldsymbol{s}_{1},k_{1},\boldsymbol{s}_{2},k_{2})$ are independent of $y_{1}$ and $y_{2}$, and $\mathcal{P}(x_{1},y_{1})$ is independent of $y_{1}$, we expect a solution in which $u_{k_{1}}(t,\boldsymbol{s}_{1})$ is independent of $y_{1}$. Integrating over $y_{2}$ gives

[TABLE]

where

[TABLE]

We solve (30) with no-flux conditions at $x=0$ and $x=1$ and initial condition

[TABLE]

because initially there are $10^{3}$ particles placed uniformly at random in $[0,1/10]\times[0,1]$.

In Fig. 4, we show a comparison between a Monte-Carlo simulation of $u_{k}(t,x)$ using Algorithm 1 and a numerical solution of equation (30) using a second-order central difference finite-volume method (FVM) in space and a fourth-order Runge–Kutta scheme in time. We use the parameter values $\mu=3/4$, $\sigma=1/4$, $\mathcal{J}=500$, and $\epsilon=0.1$.

We observe that the distribution function $u_{k}(t,x_{1})$ is nontrivial, and that the reduced model does a good job of capturing the empirical distribution that we obtain from Monte-Carlo simulations. The KS probability for the distributions in Fig. 4 is $\rho_{\text{KS}}\approx 3.2\times 10^{-2}$.

VI.3 Limit of Large Networks

We expect our mean-field assumption to be more accurate for larger networks. We now briefly investigate this hypothesis in the context of degree distributions. (Naturally, it is also relevant to consider this hypothesis for other quantities.)

For a network with $n$ nodes, $u_{k}=O(n)$, so we see that equation (25) converges to a sensible limit as $n\rightarrow\infty$ if $\mathcal{C}=O(1/n)$ and $\mathcal{D}=O(1)$. To investigate the convergence for large $n$, we consider networks in which the number of nodes is constant in time (i.e., $\mathcal{J}\equiv 0$). We first consider networks in which there is no edge deletion (i.e., $\mathcal{D}\equiv 0$). We take the edge-creation rate to be $500/n$ times that given in equation (27). We choose all other parameters as in Section VI.2.

In the top panel of Fig. 5, we show the mean degree distribution sampled over $10^{6}/n$ realizations of Algorithm 1 for $n=125$, $n=500$, $n=2000$, and $n=8000$. We also show the degree distribution calculated by solving the IPDE (30). Qualitatively, this figure supports the hypothesis that, at least in terms of degree distribution, the mean-field approximation is more accurate for larger networks.

We now introduce edge deletion and choose $\mathcal{D}$ to be given by equation (28). We increase the rate of edge creation slightly by taking the edge-creation rate to be $750/n$ times that given in equation (27). We show the resulting mean degree distribution sampled over $10^{6}/n$ realizations of Algorithm 1 for $n=125$, $n=500$, $n=2000$, and $n=8000$ in the bottom panel of Fig. 5(bottom) along with the degree distribution that we calculated by solving the IPDE (30).

VI.4 One-Step Network Creation Versus Kinetic Approximations

Dynamic models of network creation can provide an alternative to one-step network creation. For example, in the standard $G(n,p)$ Erdős–Rényi (ER) model Newman_2010 , one specifies that a network has $n$ nodes and that each pair of nodes is connected with independent, constant probability $p\in(0,1)$. This leads to a binomial degree distribution $\text{Bin}(n-1,p)$, which becomes the Poisson distribution $\text{Pois}(np)$ in the limit $n\to\infty$ with fixed $np$. Reference Krapivsky_2010 discussed an alternative, dynamic approach to the ER model in which an initially unconnected network has $n$ nodes and each node connects to other nodes uniformly at random at a specified rate. In the $n\to\infty$ limit, one can solve a master equation to obtain a Poisson degree distribution, which coincides with the standard model when halted at a specific time (depending on the edge-creation rate). We note also the work by Krioukov and Ostilli Krioukov_2013 , who showed that certain equilibrium ensembles create the same distribution of graphs as nonequilibrium ensembles.

References Boguna_2003 ; Boguna_2003_b discussed a one-step network-creation model that allows nodes to have an associated state. Nodes $i$ and $j$ have randomly distributed latent social variables $\boldsymbol{s}_{i}$ and $\boldsymbol{s}_{j}$, and these nodes are adjacent to each other with probability $r(\boldsymbol{s}_{i},\boldsymbol{s}_{j})$. References Boguna_2003 ; Boguna_2003_b then used a mean-field approximation to derive a formula for the degree distribution (depending on the probability distribution of the social variables and on the function $r$) of the network as the number of nodes tends to infinity.

One can use our kinetic approach as an alternative to one-step creation to examine such a scenario. To demonstrate this, we consider the case investigated in Boguna_2003 , where the latent social variable is a positive scalar $\boldsymbol{s}=h\in[0,h_{\max}]$ and

[TABLE]

for constants $b$ and $\alpha$. We solve equation (IV.2) with an initial condition corresponding to $n$ unconnected nodes that are distributed uniformly at random in a state space $[0,h_{\max}]$. No new nodes enter the system (i.e., $\mathcal{J}\equiv 0$), there is no edge deletion (i.e., $\mathcal{D}\equiv 0$), and the edge-creation rate is $\mathcal{C}=r(h_{i},h_{j})$. Thus, at time $T_{\text{end}}=1$, the expected number of edges between node $i$ and node $j$ is $r(h_{i},h_{j})$. In Fig. 6, we see for parameter values $n=1000$, $\alpha=3$, and $b=1/2$ that the degree distribution given by our mean-field model at $t=T_{\text{end}}$ matches very closely with the analytical formula of Boguna_2003 . We solved the IPDE using a method of lines with a spatial discretization of $\Delta h=1/4$ and a fourth-order Runge–Kutta scheme in time.

VI.5 Final Remarks: Moment

Closure

Despite the fact that equations (IV.2) and (25) are of much lower dimension than equations (7) and (23), each of the former equations is still an infinite system of IPDEs, and they may still be too expensive to solve numerically. One can make further approximations by considering moments of the density $u_{k}(t,\boldsymbol{s})$ with respect to the degree $k$:

[TABLE]

In particular, $\mathcal{M}^{(0)}$ gives the number density of nodes (for which we obtain a closed equation if $\mathcal{J}\equiv 0$ and $\mathcal{L}^{(1)}$ does not depend on the network structure), and $\mathcal{M}^{(1)}$ gives the number density multiplied by the mean degree of a node with state vector $\boldsymbol{s}$. It is only possible to obtain a closed system of equations for particular choices of $\mathcal{C}$ and $\mathcal{D}$. In general, one needs to truncate the hierarchy of moment equations and then apply another closure assumption by positing an expression for a high-order moment in terms of lower-order moments Hillen_2003 ; Hillen_2004 ; Kuehn_2015 .

VII A Model for Osteocyte Network Formation

VII.1 Osteocytes

An osteocyte KnotheTate_2004 is a dendritic cell in both cortical bone (dense, weight-bearing bone) and trabecular bone (flexible, highly vascular bone). The protrusions (i.e., dendrites) of the cell are known as processes, and they form a communication network between the osteocytes and cells on the bone surface. To avoid confusion with the word processes'', which is common jargon in physics, we henceforth refer to *processes* as dendrites''. Osteocytes are densely packed in bone. They occupy spherical spaces called lacunae, and their dendrites occupy tunnels called canaliculi Carter_2014 ; Piattelli_2014 . The dendrites enable communication via gap junctions, across which signaling molecules can diffuse. The precise purpose of this communication network is not known, but there are many conjectures. It has been suggested that the exchange of signaling molecules relates to skeletal unloading, fatigue damage, and estrogen deficiency Jilka_2013 . The range of signaling molecules that have been detected is vast, and many also arise in the regulation of other organs. These include Receptor Activator of Nuclear Factor Kappa-B Ligand (RANKL) Jilka_2013 , Vascular Endothelial Growth Factor (VEGF) Jilka_2013 , Parathyroid Hormone (PTH) Xiong_2014 , calcium ions ($\text{Ca}^{2+}$) Ishihara_2012 , and Sclerostin Sapir_Koren_2014 . Additionally, there is a thin layer of fluid around the osteocyte network. Perturbations to osteocyte-network organization can impact both fluid flow and diffusion, and they thereby allow mechanosensation and signaling Kerschnitzki_2013 .

Because of the location of osteocyte networks in bone, it is difficult to examine them experimentally. In Fig. 7, we show the outcome of applying an obtrusive experimental technique to view osteocytes after nearby mineral has been dissolved. Three-dimensional imaging data is now available by using confocal microscopy Kamioka_2001 ; Sugawara_2011 ; Sugawara_2005 , and such work has led to the identification of some structural features of osteocyte networks. Identified features include the mean number of dendrites that protrude from each osteocyte Sugawara_2005 and mean lengths of a canicular network Buenzli_2015_b . Other work has reported that high-density networks correlate positively with high bone quality. Note that the mineral matrix has an orientation (from dendrites and collagen), and bone quality is associated with the level of organization of this mineral matrix Kerschnitzki_2013 .

VII.2 Formation Process

On the bone-tissue interface, two cell types are actively involved in the bone-formation process: osteoblasts and osteoclasts. Osteoblasts form a layer on the bone surface and secrete the osteoid bone matrix. The larger multi-nucleated osteoclasts subsequently resorb the bone matrix Robling_2006 . Osteoblasts also express RANKL and osteoprotegerin (OPG), which promote and inhibit the bone resorption by osteoclasts, respectively. This is one example for how osteoblasts tightly regulate bone formation and destruction. As the osteoblasts produce the calcium matrix, occasionally they become embedded within the bone. These osteoblasts then change morphology to become star-shaped osteocytes.

Osteoblasts originate from mesenchymal cells and have one of four possible fates: undergo apoptosis (approximately 65%), become embedded in bone as osteocytes (approximately 30%), transform into inactive osteoblasts and become bone-lining cells, or transdifferentiate into cells that deposit chondroid bone Franz-Odendaal_2006 . Upon some signaling event, osteoclasts arrive at a bone and the osteoblasts move aside. The osteoclasts then burrow into the bone; as they do so, they resorb some of the osteocyte matrix. It has been suggested that after an osteocyte undergoes apoptosis, pro-osteoclastogenic signals are released by the osteocyte's neighbors in the network Kennedy_2014 . A trail of osteoblasts then follows the osteoclasts and secretes new bone matrix 333The unit that consists of osteoblasts following osteoclasts is known as a “Bone Multicellular Unit” (BMU)., although some of these get left behind to become osteocytes. Kamioka et al. suggested that osteoblasts are incorporated into a network by osteocytes extending their dendrites towards the osteoblast layer Kamioka_2001 .

Thus far, we have discussed three types of bone cells: osteoblasts, osteoclasts, and osteocytes. For at least the osteoblast-to-osteocyte cell transition, biologists have subdivided the process of cell differentiation to include eight phenotypes: (i) preosteoblast; (ii) preosteoblastic osteoblast; (iii) osteoblast; (iv) osteoblastic osteocyte; (v) osteoid-osteocyte (i.e., Type-II preosteocyte); (vi) Type-III preosteocyte; (vii) young osteocyte; and (viii) old osteocyte Franz-Odendaal_2006 . Additionally, the secretion of bone occurs as two steps: first osteoid is deposited as a scaffold, and then mineralization occurs to confer strength. Stages (iv)–(vi) are cells after the deposition front but before the mineralization front; they are surrounded by a non-mineralized osteoid matrix. (In other words, there is scaffold around them.) Stages (vii)–(viii) are cells whose volume has depleted (by reduction in the endoplasmic reticulum and Golgi apparatus) and are in mineralized bone. The diagram in Fig. 8 shows the bone-formation step. Although it is potentially useful to consider all of the above phases (defined from osteogenic markers), we are interested only in the structure of a mature osteocyte network [stages (vi)–(viii)], so we will make drastic simplifications.

VII.2.1 Bone Metastasis

Advanced prostate, breast, and lung cancer can metastasize to bone Roudier_2003 ; Zhang_2013 . In pathological bone, the highly regulated bone-remodeling signaling pathway is disrupted. A particularly painful symptom is net bone formation in some regions and simultaneous weakening in other areas Pivonka_2010 . Small bone lesions can also develop.

In a cancerous microenvironment, transforming growth factor beta (TGF$\beta$) expressed by tumor cells promotes excessive osteoblast growth Lee_2011 . The overexpression of TGF$\beta$ is only one of the many ways in which tumors can interfere with bone formation; TGF$\beta$ targeted cancer therapy has only been marginally successful Connolly_2012 .

Ordinary differential equation (ODE), and hence non-spatial, models for cell populations in healthy bone that incorporate osteoblast, osteoclast, and osteocyte populations were developed in Graham_2013 ; Moroz_2006 ; Wimpenny_2007 . Partial differential equation (PDE) models of healthy bone remodeling include Ryser_2010 ; Ryser_2009 , and these models were adapted subsequently for cancerous bone in Ryser_2012 . Mechanically-focused models that capture stresses and strains in bone have also been explored Pivonka_2010 ; Rejniak_2011 ; Oers_2014 . A few of these models consider osteocyte density, but none of them explore network structure.

A seemingly unexplored area is the investigation of osteocyte network morphology in the presence of cancer. There is evidence in Stinson_1975 ; Eisenberger_2008 that for myeloma and (benign) osteoma, osteocytes are exceptionally spherical and have shorter, distorted dendrites that are reduced in number. A contrasting osteocyte network with unregulated excessive dendritic growth (and hence larger numbers of dendrites) was observed in the presence of osteogenic sarcoma Stinson_1975 .

In Section VII.3, we develop a simple model that incorporates network properties into the bone-formation process. For some cancer types, it is known in part how bone formation is affected (e.g., there is increased osteoblast proliferation). Using such a model, it may be possible to connect a change in the bone-formation process to the properties of the resulting osteocyte network. Qualitatively, one would then be able to suggest which osteocyte network phenotype is promoted by a particular cancer (e.g., stunted dendrite growth, excessive dendrite growth, etc.). Conversely it may be possible to infer changes in the bone-formation process from observations of osteocyte network structure.

VII.3 Model of Osteocyte Network Growth

The model that we develop in this section builds on the work of Buenzli Buenzli_2015 . We avoid modeling the full complications of the biology (e.g., cell sub-classifications, proteins, etc.) and consider only osteoblasts and osteocytes. By examining a simple model, we hope to gain insight into how osteocyte network structure may depend on measurable quantities, while preserving a minimalist approach.

The osteocyte network occupies an expanding domain $\Omega(t)\subset\mathds{R}^{d}$ with a boundary $\partial\Omega(t)$ that moves with normal velocity $v(t)$. The nodes in the network are osteocytes, which each have an associated position in $\mathds{R}^{d}$. The undirected edges are the dendrites between them. We allow multiedges, which correspond to multiple dendritic connections between the same pair of osteocytes, and they are observed in practice Kerschnitzki_2013 . We suppose that the positions of the osteocytes are time-independent, so $\mathcal{L}^{(1)}\equiv 0$.

Motivated by the observation that osteoblasts differentiate to osteocytes near the mineral front Franz-Odendaal_2006 — where there is less mineralization — we suppose that the rate at which an osteocyte creates connections to others is governed by the local bone mineral density $m(t,\boldsymbol{x})\in\mathds{R}^{+}$.

Our model consists of two processes: (1) domain expansion and (2) edge creation within semi-mineralized bone. Buenzli's model Buenzli_2015 consists entirely of domain expansion, whereas we also incorporate a network structure.

VII.3.1 Domain Expansion

We suppose that osteoblasts are encased within the mineral matrix and become osteocytes with degree $k=0$ at a rate of $D_{\text{burial}}(t,\boldsymbol{x})\rho_{\text{Ob}}(t,\boldsymbol{x})$ for $\boldsymbol{x}\in\partial\Omega(t)$, where $D_{\text{burial}}(t,\boldsymbol{x})$ is the probability per unit time of an osteoblast joining the matrix and $\rho_{\text{Ob}}(t,\boldsymbol{x})$ is the surface density of osteoblasts. (We take $\rho_{\text{Ob}}(t,\boldsymbol{x})$ to be given.) Equation (IV.2) is then

[TABLE]

where

[TABLE]

Following Buenzli_2015 , we suppose that the (outward) normal velocity of the interface is

[TABLE]

where $\kappa_{\text{form}}(t,\boldsymbol{x})$ is the volumetric rate at which osteoblasts form the mineral matrix.

VII.3.2 Edge Creation in Semi-Mineralized Bone

We suppose that the mineral density $m(t,\boldsymbol{x})$ of the matrix can vary in both space and time, as mineral is produced by both osteoblasts on the deposition front and osteocytes behind it. We suppose that there is a maximum mineral density $C_{m}$ that corresponds to fully mineralized bone, that osteocytes produce mineral at a rate of $r_{\text{exc}}^{(\text{Cy})}(1-{m}/{C_{m}})$ per cell, and that osteoblasts produce mineral at a rate of $r_{\text{exc}}^{(\text{Ob})}$ per cell, where $r_{\text{exc}}^{(\text{Cy})}$ and $r_{\text{exc}}^{(\text{Ob})}$ are constants. Therefore,

[TABLE]

where $f(t,\boldsymbol{x})$ [given by equation (1)] is the density of osteocytes. Note that $\kappa_{\text{form}}$ is the volume of matrix produced per osteoblast per unit time, whereas $r_{\text{exc}}^{(\text{Ob})}$ is the mass of matrix produced per osteoblast per unit time, so the density of matrix formed by the osteoblasts is $r_{\text{exc}}^{(\text{Ob})}/\kappa_{\text{form}}$, which should be no larger than the maximum mineral density $C_{m}$. Note additionally that because $D_{\text{burial}}$ is the rate at which osteoblasts join the matrix and $\kappa_{\text{form}}$ is the rate at which matrix volume is produced, the ratio $D_{\text{burial}}/\kappa_{\text{form}}$ is the density of newly formed osteocytes.

We now model how edges (i.e., dendrites) form. We suppose that the rate of edge creation depends on the mineral density $m$. In Kamioka_2001 , it was suggested that osteocytes grow dendrites towards the osteoblast layer. This suggests that dendrites grow in the part of the domain that is not fully mineralized. To construct a simple model in which edges are less likely to form as mineral density becomes larger, we let the rate of edge creation between nodes at positions $\boldsymbol{x}$ and $\boldsymbol{y}$ be

[TABLE]

where $g$ is monotonically decreasing and vanishes at infinity, encoding the fact that short edges are much more likely to form than long edges. Equation (38) takes into account the mineral density only at the two endpoints $\boldsymbol{x}$ and $\boldsymbol{y}$. One can formulate more complicated models in which $\mathcal{C}$ depends on (for example) a line integral of $m$ between $\boldsymbol{x}$ and $\boldsymbol{y}$. However, given the simplifications and modeling assumptions that we have already made, we do not consider such complicated edge-creation models.

With our model assumptions, equation (IV.2) becomes

[TABLE]

Equation (39) is coupled through $\mathcal{C}$ to equation (37) for the mineral density, and the domain $\Omega(t)$ evolves according to equation (36).

Summing (39) over $k$ yields

[TABLE]

which is identical to equation (8) in Buenzli_2015 . Note also that equation (39) does not include feedback between network structure and the burial rate $D_{\text{burial}}$, though one can extend the model to incorporate such coupling.

VII.4 Traveling-Wave Solution

If $D_{\text{burial}}$, $\kappa_{\text{form}}$, and $\rho_{\text{Ob}}$ are constant, then equations (36), (37), and (39) admit a traveling-wave solution that corresponds to the sustained creation of new bone. Such a solution may give an indication of the local behavior near growing bone, and we can use it to determine how the properties of the bone depend on the parameters in the model.

For a one-dimensional traveling wave that moves with a constant speed of $v=\kappa_{\text{form}}\rho_{\text{Ob}}$, we can (without loss of generality) take the domain to be $\Omega(t)=(-\infty,vt)$. We transform to coordinates that move with the wave by writing $z=x-vt$, and we seek a solution in which $u_{k}$ and $m$ depend only on $z$. Equation (39) becomes

[TABLE]

Summing over $k$ gives

[TABLE]

so

[TABLE]

where $H(x)$ is the Heaviside function. Equation (42) is identical to equation (13) in Buenzli_2015 when $D_{\text{burial}}$ and $\kappa_{\text{form}}$ are constant. Using (42) in equation (41) gives

[TABLE]

where

[TABLE]

Equivalently, one can write equation (43) as

[TABLE]

whose solution is

[TABLE]

The density of newly-formed osteocytes $b$ sets the scale for $u_{k}$.

Before we can solve for $\lambda$, we first need to solve for $m(z)$. Solving equation (37) yields

[TABLE]

We see that the mineral density varies from that produced by osteoblasts ${r_{\text{exc}}^{(\text{Ob})}}/{\kappa_{\text{form}}}$ to the maximum mineralization $C_{m}$ over a length scale of

[TABLE]

We illustrate this behavior in Fig. 9(top).

Before we can evaluate $\lambda$, we need to choose a form for $g$. Suppose first that $g$ is a constant, and let's write $g\equiv\beta$. In principle, this may allow long edges to form — recall (see equation (38)) that we argued that $g$ should decay at infinity to ensure a low probability for long edges to form — but we see from equation (46) that even with $g$ identically constant, $C(z,z^{\prime})$ approaches [math] exponentially fast over the length scale $L$, which therefore sets the scale for the maximum edge length.

Inserting equation (46) in equation (45) gives $\lambda=\lambda_{\infty}\left(1-\mathrm{e}^{z/L}\right)$, where

[TABLE]

which indicates that the mean degree varies over the same length scale from a value of [math] for newly-formed osteocytes to a value of $\lambda_{\infty}$ deep within a bone.

In the bottom panel of Fig. 9, we illustrate the traveling-wave profile of the degree distribution by plotting $(1/b)\sum_{k=0}^{K}u_{k}(z)$ for several values of $K$. The differences between these curves indicate the proportion of osteocytes of each degree. As expected, there is a region at the front of the wave in which the mean degree of the osteocytes is lower, while the LSDD approaches a stationary distribution far behind the front.

If, instead of taking $g$ to be constant, we instead choose $g(z,z^{\prime})=\beta\mathrm{e}^{-|z-z^{\prime}|/l}$ or $g(z,z^{\prime})=\beta\mathrm{e}^{-|z-z^{\prime}|^{2}/l^{2}}$, we obtain qualitatively similar results, provided $l>L$. For example, in the first case,

[TABLE]

For $l\gg L$, we see that $\lambda\to\lambda_{\infty}\left(1-\mathrm{e}^{z/L}\right)$. For $l\ll L$, we see that $\lambda\to\lambda_{\infty}(l/L)\left(1-\mathrm{e}^{2z/L}\right)$, so the length scale over which degree varies is halved and the mean degree smaller by the factor $l/L$.

VII.5 Parametrization and Interpretation

One can infer representative values of some of the parameters and variables in our model from existing experimental data. For example, in the review paper  Buenzli_2015_b , the authors calculated the number of osteocytes per $\text{mm}^{3}$ to lie in the range $19000$–$28500\,\text{mm}^{-3}$; the surface density $\rho_{\text{Ob}}$ of osteoblasts per $\text{mm}^{2}$ was calculated in Buenzli_2014 (using data from Ref. Marotti_1975 ) to be in the range $2000$–$10000\,\text{mm}^{-2}$; and values of $D_{\text{burial}}$ and $\kappa_{\text{form}}$ were given in Buenzli_2015 (using data from Refs. Marotti_1975 ; Hannah_2010 ) for bone near a Haversian canal (which occurs only in cortical bone). The number of dendrites that protrude from an osteocyte yielded a mean degree of $\langle k\rangle\approx 52.7$ in Sugawara_2005 . As suggested in Kamioka_2001 ; Sugawara_2011 ; Kerschnitzki_2013 , we expect the number of unique neighbors of an osteocyte to be smaller by an order of magnitude. Measurements of the speed of the mineralization front of bone creation were given in Araujo_2014 . Numerous other parameters (e.g., the rate $r_{\text{exc}}^{(\text{Cy})}$ at which osteocytes produce mineral and the maximum mineral density $C_{m}$) are unknown, and the values above were reported in only a few papers, so it is not yet possible to make quantitative predictions with our model.

Despite the dearth of knowledge about parameter values, we can use our model to make qualitative predictions of the effect of varying each parameter. The review Logothetis_2005 summarized how different cancer types interfere with healthy bone remodeling. In particular, the review details how factors produced by prostate cancer cells lead to net bone formation due to increased levels of proliferation and differentiation of osteoblasts. In our model, an increased level of differentiation corresponds to an increased burial rate $D_{\text{burial}}$. Increasing $D_{\text{burial}}$ in the model increases the number density of osteocytes $f$, but it decreases the mean degree $\lambda$. In other words, it decreases the mean number of dendrites that leave a cell body. In the model, this occurs because an increase in osteocyte density leads to an increase in mineralization, which makes dendrite formation less likely.

It was shown in Stinson_1975 ; Eisenberger_2008 that for myeloma and benign osteoma—two other cancers that can cause net bone formation Bataille_1990 —osteocytes are rather spherical with shorter, distorted dendrites that are fewer in number in comparison to those in healthy tissue. Although these articles do not comment on osteocytes density or on mineralization, and it is crucial to be careful to avoid over-interpreting such limited data, our model does suggest one possible mechanism for this change in morphology.

VIII Conclusions

We introduced a model for evolving spatial networks, and we used a mean-field approximation to reduce the dimension of its governing hierarchal Fokker–Planck equations. Specifically, by defining a local state degree distribution, we derived IPDEs (IV.2) and (25) to describe an evolving spatial network that includes evolution of the position of nodes (or some more general state vector), edge creation, edge deletion, and new node creation that occur at prescribed rates. Our approach generalizes commonly-studied master-equation approaches by including a state space so that we can examine spatial networks.

To illustrate the potential utility of our IPDEs in applications, we examined growing osteocyte networks in bone. Although we employed a very simplistic model, we were able to use it predict relationships between biological parameters and network structure. Our approach provides a starting point for examining spatial networks in biology and other fields. In the future, we hope to compare predictions to experimental data after making a model more faithful to the biology.

IX Acknowledgements

J. P. T.-K. received funding from the EPSRC under grant reference number EP/G037280/1. We thank Pascal Buenzsli, Mikko Kivelä, Andrew Krause, SeWook Oh, Jan Poleszczuk, and particularly Paul J. Dellar for helpful discussions.

Appendix A Candidate Algorithm for the Kinetic Network

Model

In Algorithm 1, we give pseudocode for our simulations of evolving spatial networks. In the main text, we summarized our model in Table 1.

When simulating Algorithm 1, we use a small time step $\Delta t$, so $0<\Delta t\ll 1$. We also specify the following ordering of events: edge creation; edge deletion; state update; and then new particles are allowed to enter the system. This specification is arbitrary, and obviously it is desirable that any reordering of these events becomes inconsequential as $\Delta t\rightarrow 0$. For our numerical experimentation using the examples in Section VI, this indeed appears to be the case. For the simulations that we reported in Section VI, we chose a time step $\Delta t$ to be sufficiently small that a reordering of events has no discernible impact.

The update rule for event (iv) in Table 1—namely, the evolution of the states of the individual nodes—depends on the particular process that we consider. We write $\boldsymbol{s}_{i}(t+\Delta t)=\mathscr{D}(\boldsymbol{s}_{1}(t),\dots,\boldsymbol{s}_{N(t)}(t),\Delta t)$, where $\mathscr{D}$ arises from the time-discretization of the underlying process. For example, one can use an Euler–Maruyama method or the Milstein method for the SDE (4) Higham_2001 ; and one use Störmer–Verlet schemes for the ODEs in equation (3).

One can devise a efficient simulation algorithm for situations in which edge creation or deletion do not depend on the state of nodes. In that case, one can use an event-driven algorithm, such as a Gillespie algorithm Gillespie_1977 , for creation and deletion events.

Appendix B Additional Details of the Model

Derivation: No Edge Deletion

In this appendix, we fill in the details of the derivation in Section IV.2. When we sum over the degrees and integrate over the states of particles 2 through $n$, multiply by $n$, and sum over $n$, the first term on the RHS of equation (7) gives

[TABLE]

For $i>1$, each individual term appears once in the positive sum and once in the negative sum; they thus cancel each other out. The remaining terms are

[TABLE]

Because $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$ is invariant with respect to particle relabeling, we can relabel $(\boldsymbol{s}_{j},k_{j})\leftrightarrow(\boldsymbol{s}_{2},k_{2})$ in each term in the sum over $j$ to obtain

[TABLE]

where the last line follows from equation (14). Consequently, we can write equation (48) as

[TABLE]

For the remaining terms on the RHS of equation (7), we sum over the degrees and integrate over the states of particles 2 through $n$, multiply by $n$, and sum over $n$ to obtain

[TABLE]

For the term $i=1$, we are summing and integrating over all arguments of $\mathcal{F}_{n-1}^{\,\vec{k}_{n}^{i-}}(t,\vec{s}_{n}^{i-})$, so

[TABLE]

where the last line follows from equation (5). For each term $i>1$, we use the invariance of $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$ with respect to particle relabeling to swap particle $i$ with particle $n$ to give

[TABLE]

which cancels with the remaining term in equation (51).

Appendix C Additional Details of the Model Derivation: Edge

Deletion

In kinetic theory, the first question to address when deriving a reduced model is which variables to retain in the model and which to integrate over. In Section IV, each variable was associated with a node in a network, and it was natural to integrate over all nodes but the first. We thus retained the state and degree of node 1 as independent variables. We could have reduced the model further by subsequently integrating over either the state or degree of node 1.

When considering edge deletion, it is much more difficult to associate the independent variables with individual nodes, and each entry of an adjacency matrix is associated with a pair of nodes. Consequently, it is not obvious which variables are natural to retain in a reduced model and which variables should be integrated out. To facilitate a direct comparison of the reduced model including edge deletion with the reduced model of Section IV, we again retain the state and degree of node 1 as independent variables. We thus sum over all entries of the adjacency matrix for which node 1 has degree $k_{1}$. Additionally, as before, we integrate over $\boldsymbol{s}_{n}^{(2)}$, multiply by $n$, and sum over $n$.

Because the operators $\mathcal{L}^{(n)}$ that we are considering do not depend on network structure, the approximation of the LHS of equation (23) proceeds as in Section IV.2. For the edge-creation term on the RHS of equation (23) we find, as in Section B, that for $i>1$, each term appears once in the positive sum and once in the negative sum; these terms thus cancel. For the remaining terms (for which $i=1$), we exploit the invariance of $\mathscr{F}_{n}^{\,A_{n}}(t,\vec{s}_{n})$ with respect to particle relabeling. Specifically, we relabel $j\leftrightarrow 2$ in each term in the sum over $j$ (i.e., swapping rows and columns of $A_{n}$) to obtain

[TABLE]

where the Kronecker $\delta$ enforces the degree condition. The 2-particle LSDD is

[TABLE]

so one can write equation (52) as equation (50), and the analysis proceeds as in Section B.

Let's now consider the edge-deletion terms. As with the edge-creation terms, for $i>1$, each term appears once in the positive sum and once in the negative sum; these terms thus cancel each other. For the remaining terms (for which $i=1$), relabeling $j\leftrightarrow 2$ in each term in the sum over $j$ yields

[TABLE]

where

[TABLE]

There is now a new closure problem, as we need to relate $U_{k_{1},k_{2}}^{(2)}(t,\boldsymbol{s}_{1},\boldsymbol{s}_{2})$ to known variables. We write

[TABLE]

where $\alpha$ is the expected number of edges between nodes 1 and 2, given that these nodes have degrees $k_{1}$ and $k_{2}$, respectively. If there are $m$ edges in total, there are $2m$ stubs, of which $k_{1}$ are at node 1 and $k_{2}$ at node 2. One can approximate the probability that a given edge connects nodes 1 and 2 as

[TABLE]

so a reasonable closure assumption (reminiscent of a configuration model and hence with similar associated assumptions Fosdick2016 ) for the expected number of edges between nodes 1 and 2 is

[TABLE]

where $E[N]$ is the expected number of nodes and $\langle k\rangle$ is the mean degree. Using a mean-field approximation for $u_{k_{1},k_{2}}^{(2)}(t,\boldsymbol{s}_{1},\boldsymbol{s}_{2})$, we can then close the edge deletion term by writing

[TABLE]

Finally, we consider the node-creation term. The term $i=1$ gives

[TABLE]

where the last equality follows by equation (20). For each term with $i>1$, we can use the invariance of $\mathcal{F}_{n}^{\,\vec{k}_{n}}(t,\vec{s}_{n})$ with respect to particle relabeling to swap particle $i$ with particle $n$. As in our prior calculations, we then find that all of these terms cancel each other out.