Random walk models for the propagation of signalling molecules in one-dimensional spatial networks and their continuum limit

The propagation of signalling molecules within cellular networks is affected not only by network topology, but also by the spatial arrangement of cells in the networks. Understanding the collective reaction–diffusion behaviour in space of signals propagating through cellular networks is an important consideration, for example, for regenerative signals that convey positional information. In this work, we consider stochastic and deterministic versions of random walk models of signalling molecules propagating and reacting within one-dimensional spatial networks with arbitrary node placement and connectivity. By taking a continuum limit of the random walk models, we derive an inhomogeneous reaction–diffusion–advection equation, where diffusivity and advective velocity depend on local node density and connectivity within the network. Our results show that large spatial variations of molecule concentrations can be induced by heterogeneous node distributions. Furthermore, we find that noise within the stochastic random walk model is directly influenced by node density. We apply our models to consider signal propagation within the osteocyte network of bone, where signals propagating to the bone surface regulate bone formation and resorption processes. We investigate signal-to-noise ratios for different damage detection scenarios and show that the location of perturbations to the network can be detected by signals received at the network boundaries.


Introduction
Biological tissues are composed of complex arrangements of cells in space that enable them to perform a variety of functions.To understand how these functions emerge from cell-cell interactions taking place in these complex arrangements, new quantitative approaches based on networks have been developed [1,2,3,4].These networks represent cell-cell communication pathways due for example to the adjacency of cells in plants and confluent epithelial cell layers [5,6,7], or due to branching cellular dendrites connecting near and far cells, such as in the neural network in the brain [8,9] and in the osteocyte network in bone [10,11,12,13].Distributed networks also arise in cardiovascular, respiratory, plant vascular and root systems [14].In all these biological systems, understanding how network architecture influences transport properties of signals propagating through these networks is important to assess tissue functional behaviour [15].For example, transport properties may affect the efficiency of nutrient delivery and waste product removal in surrounding tissues.In plants, signalling through the cellto-cell communication network regulates growth, patterning, and defense mechanisms [16,7].
Random walks are an effective way to describe stochastic reaction-diffusion processes occurring on networks [17,18].They have been used to describe disease spread in social networks, the spread of opinions in voting models [19], and transportation systems [20,21].Several common reaction-diffusion continuum models, such as Fisher-KPP and activatorinhibitor models, are generalised to networks using discrete Laplacians to investigate signal propagation [22,23,24,25] and spatio-temporal patterns induced by dynamic instabilities [26,27,28,29,30].
Most studies of signal propagation in networks focus on the impact of the topological graph structure.However, many biological networks are embedded in physical space [31].The function of some of these networks relies on signals transmitting positional information, which depends not only on network topology, but also on the physical location of the network's nodes and edges.The osteocyte network in bone, for example, senses local mechanical deformation of bone tissue and detects micro-damage to direct bone adaptation and bone regeneration where needed [32,33,34].This occurs through signalling molecules propagating through the osteocyte network, or through the lacuno-canalicular network in which osteocytes reside, from within bone tissue to the bone surface.Recent mathematical models have simulated fluid flow in lacunocanalicular networks of mouse tibia using circuit theory on weighted graphs [35,36,37].While these simulations take into account precise network information from micro-CT scans, they provide limited mathematical understanding of the effect of network architecture on average signal transmission at the tissue scale.They also do not take into account the effect of chemical reactions such as creation and degeneration of signalling molecules and the importance of stochasticity in these processes.
In this paper, we investigate random walk models of signal propagation in one-dimensional spatial networks where node placement and connections between nodes is arbitrary (Section 2.1).We assess the collective reaction-diffusion behaviour of these models by deriving a continuum limit in which the physical spacing between the nodes is small (Section 2.2).The theory of random walks on regular lattices is well-developed and known to be described by reactiondiffusion phenomena [38,39,40,24].Here, we show that on irregular networks, the continuum limit leads to an inhomogenous reaction-diffusion-advection partial differential equation for molecule concentration.The diffusivity and advective velocity in this equation depend explicitly on local node density and connectivity.We also formulate an equivalent deterministic positionjump model that captures the average behaviour of the stochastic random walk model and provides fast and accurate solutions to the continuum equation [41] (Section 2.1).To explore the role of the spatial arrangement of nodes, we perform numerical simulations of random walks on a series of irregular spatial networks and compare these results with numerical solutions of the reaction-diffusion-advection equation (Section 3).These comparisons emphasise the importance of the spatial structure of networks for quantities related to signal propagation, such as transport and signal-to-noise ratios.Finally, we provide an application of these models to the osteocyte network in bone (Section 3.3).The reaction-diffusion-advection equation is of particular interest in this system to capture coarse-grained behaviour of signal propagation through this dense, large spatial network [12].

Methods
We first present the stochastic and deterministic variants of the random walk model in Section 2.1.The evolution of average node occupancy is governed by Eq. (5), and a relationship between the variance of the molecular concentration and the node density is established in Eq. (13).We then derive a continuum limit of the deterministic model in Section 2.2, which leads to a reaction-diffusion-advection equation (Eq.( 31)) with a diffusivity and an advective velocity that depend explicitly on the network's node density and connectivity (Eqs ( 33)-( 34)).

Stochastic and deterministic discrete models
We consider a general discrete-time random walk model for the propagation of signalling molecules in a one-dimensional spatial network.The network's nodes are distributed arbitrarily along the x -axis at positions x i , i = 1, . . ., K , that are ordered by increasing value (Fig. 1a).Each node i is assumed to hold a number of signalling molecules N i (t ) at time t .We define the territory of node i , ∆x i , as half the distance to node i − 1 plus half the distance to node i + 1: This spatial network has a node density at x i given by The concentration n i (t ) of molecules at x i at time t is defined as Each node i is assumed to be connected to an arbitrary number of nodes of the network among all available nodes j = 1, . . ., K .We introduce q s i as the probability for a signalling molecule to make a jump of size s from node i during a time increment ∆t , where s can be positive, negative, or zero, i.e., q s i is the probability for the signalling molecule to jump from node i to node j = i + s during ∆t .The network is assumed to be directed, such that the probability of jumping from node i to node j may not necessarily equal the probability of jumping from node j to node i .
Signalling molecules may also undergo chemical reactions during the same time increment ∆t .We assume that they may be created at node i with a probability per unit time C (N i ), and eliminated at node i with a probability per unit time E (N i ).These rates depend on the number of molecules present at the node and may also vary explicitly with position and time.For ease of notation we omit these additional function arguments.The net chemical reaction rate at node i is The evolution of the stochastic model is determined by formulating the balance of molecules at node i from time t to time t + ∆t .The number of molecules at node i at time t + ∆t either comes from molecules making a jump of size s from node i − s , or is due to creation or elimination of molecules according to the chemical reaction rate.Jump size s is an offset of node indices and may be a positive or negative integer, or zero if molecules remain at node i .On average, the number of molecules at node i is expected to evolve according to where the sum is taken over all possible jump sizes s leading to node i , which depends on the connectivity of node i (Fig. 1a) [39,40].Alternatively, one may consider the network to be complete (every node is connected to every node), but jump probabilities are zero where nodes are not connected.Jumps conserve molecules, so that the sum of jump probabilities for molecules leaving internal nodes must satisfy: At the boundaries of the network we assume absorbing boundary conditions, so that any incoming molecules at nodes i = 1 and i = K are removed from the system, resulting in In practice, we keep a record of the number of molecules N − (t ) and N + (t ) reaching the boundary nodes i = 1 and i = K , and calculate the flux of signalling molecules leaving the network at the left and right boundaries as respectively.Both of these fluxes are non-negative.
The stochastic model presented above provides evolution rules for the stochastic variable representing the number of signalling molecules at node i at time t (see Numerical simulations below).The stochastic model is associated with a master equation [20,38,40,39,17] that describes the evolution of the probability P i (N , t ) of having N molecules at node i at time t .In this work, we are particularly interested in describing the average behaviour of signal propagation in space in the continuum limit, which is related to 〈N stoch i (t )〉 calculated as the first moment of P i (N , t ).A closed-form evolution law for the average occupancy of molecules 〈N stoch i (t )〉 can only be found in general when reaction rates are linear or under a mean-field approximation [40,42,43].Instead of using the master equation to derive an evolution law for the average, here we propose a deterministic, discrete version of the random walk model and consider its continuum limit.Both the deterministic discrete model and its continuum limit capture the average behaviour of the stochastic random walk model.The deterministic variant of the stochastic model is obtained by considering q s i to represent the fraction of molecules making a jump of size s from node i .In this case, Eq. ( 5) provides an evolution rule for a (non-integer) variable . We refer to this deterministic variant of the model as the deterministic model.

Noise
The distribution of molecules N stoch i (t ) across the network and corresponding concentration n stoch i (t ) from Eq. (3) are random variables, whose variance quantifies the expected variability of signals propagating through the network.When there are no reactions (C (N ) = 0, E (N ) = 0), molecules propagate through the network without interacting with each other.The joint probability P i 1 ,...,i N tot (t ) that at time t molecule 1 is at node i 1 , molecule 2 is at node i 2 , etc., factorises into P i 1 ,...,i N tot (t ) = P i 1 (t )P i 2 (t ) • • • P i N tot (t ), where N tot is the total number of molecules in the network and P i (t ) is the single-molecule probability of being at node i at time t .The probability distribution of N stoch i (t ) is the probability P i (N , t ) of encountering N molecules at node i at time t , given by where the sum in Eq. ( 9) is taken over all N tot !permutations σ of node indices i 1 , i 2 , . . ., i N tot .Eq. (9) shows that at each time t , N stoch i (t ) follows a binomial distribution with parameters N tot and P i (t ).The mean and variance of N stoch i (t ) are therefore given by [44] N stoch i (t ) = N tot P i (t ), (10) Var The variance of the concentration of molecules n stoch i (t ) is found by using Eq. ( 3) and (11).We first relate P i (t ) to the average concentration of molecules along x by Eqs (10) and (3), Substituting this expression in Eq. ( 11), and since we expect that

Numerical simulations
Numerical simulations of the deterministic model are performed by updating the number of signalling molecules N disc i (t ) by stepping forward in time according to Eqs (5), (7) from an initial condition N disc i (0) for i = 1, . . ., K .
Numerical simulations of the stochastic model update the number of signalling molecules N stoch i (t ) algorithmically in two steps: a transport step, and a reaction step.These steps can be considered to occur sequentially during the same time increment ∆t under the assumption that jump probabilities are independent of the number of molecules, and consequently, independent of chemical reactions that may occur during ∆t [40].The transport step introduces an intermediate quantity, N stoch i , that represents the number of molecules at node i after transport but before any reaction.This quantity is first set to N stoch i = N stoch i (t ) at all nodes i = 1, . . ., K , and updated by accounting for jumps made by the molecules during the time increment ∆t at all internal nodes i = 2, . . ., K − 1.For each molecule at node i , i = 2, . . ., K − 1, we select a target node i + s with probability q s i , decrement N stoch i by one, and increment N stoch i +s by one.Any jump that would take a signalling molecule out of the network (i + s ≥ K or i + s ≤ 1) results in the molecule being removed from the system and contributes to the flux at the boundary.
After this transport step, N stoch i is used to process the reaction step.The probability of undergoing a creation or elimination reaction during ∆t at node i is [40] We assume that ∆t is small enough that P reac (i , t ) < 1, and that at most a single reaction occurs at each node within the time increment ∆t .If a reaction at node i occurs, then the probability that it is a creation event is in which case Conversely, the probability of the reaction being an elimination event is in which case

Continuum limit
When spacings between nodes of the network and the time increment ∆t are small, the behaviours of the stochastic and deterministic discrete models can be approximated by a continuum conservation equation for the concentration of signalling molecules.On regular lattices with even spacings and jumps to neighbouring lattice sites, the derivation of continuum limits of random walks is well established and used as a model of Brownian motion and molecular diffusion [38,17,40,39].Here, we present a derivation of the continuum limit for arbitrary one-dimensional spatial networks in which node positions and connections are arbitrary, based on the evolution equation of the deterministic model given by Eq. ( 5).
To formally take the limit as the spacing between the nodes along the x -axis approaches zero while retaining irregular spacings between nodes, we first replace the discrete mapping i −→ x i by a bĳective continuous mapping ξ → x (ξ) such that in which ξ i = i ∆ξ are evenly spaced coordinates along an axis ξ (Fig. 1b).The quantity ∆ξ represents a spatial scaling factor, so that the continuum limit is then taken by letting ∆ξ → 0 and ∆t → 0. Since node positions x i are ordered, we can assume that x (ξ) is a monotone increasing function.With the mapping x (ξ), the territory of node i is approximated by where is the metric associated with x (ξ).We define the continuous local density of nodes at position x by where ξ(x ) is the inverse function of x (ξ), such that from Eqs ( 2) and (18), ρ(x i ) ≈ ρ i .The distribution of molecules along the x -axis governed by Eq. ( 5) is mapped by the inverse function ξ(x ) into a corresponding distribution of molecules along the ξ-axis (Fig. 1).We first derive the continuum limit of the evolution equation for the concentration of molecules along the ξ-axis, where nodes are evenly spaced, and then map this evolution back onto the x -axis via x (ξ).We introduce continuous functions η(ξ, t ) and n(x , t ) representing the concentration of molecules along the ξ-axis and along the x -axis, respectively, such that The molecule concentrations along the x -axis and along the ξ-axis are related by We also introduce continuous functions of space q (ξ, s ) representing the probability for a molecule at position ξ to make a jump of size s ∆ξ along the ξ-axis, such that q (ξ i , s ) = q s i .With these definitions, Eq. ( 5) divided by ∆ξ becomes Expanding Eq. ( 23) about (ξ i , t ) to first order in ∆t and second order in ∆ξ using Assuming now that the equation holds for any value of ξ, using Eq. ( 6), and noting that are the first and second moments of the jump probabilities, respectively, Eq. ( 24) becomes where we neglected terms of order (∆t ) and (∆ξ 3 /∆t ) because they are subdominant under the assumption that the continuum limit ∆t → 0 and ∆ξ → 0 is taken such that the ratio ∆ξ 2 /∆t is constant.We also assume that jump probability biases are proportional to ∆ξ so that 〈s (ξ)〉 ∆ξ/∆t = ∆ξ 2 /∆t and both terms in the square bracket are of order (1) in the continuum limit.The above expansion in ∆ξ is valid provided node degrees are finite, or q (ξ, s ) → 0 as s → ±∞ uniformly in ξ and sufficiently fast that the moments in Eq. ( 25) are well-defined.Finally, we define and expand t ) in Eq. ( 26) to obtain Equation ( 28) is a reaction-diffusion-advection equation for the concentration of molecules along the ξ-axis, with inhomogenous diffusion D (ξ) and inhomogenous advective velocity v (ξ) that depend on node connectivity and jump probabilities via Eqs ( 27) and (25).We note here that this developement is similar to a Kramers-Moyal expansion of the master equation that leads to the Fokker-Planck or Smoluchowski equation, a partial differential equation (PDE) that governs the time evolution of the probability density function of molecule distribution in space [45,40,39,38].The advantage of not considering the master equation is that we can derive a PDE governing the average molecule concentration in the continuum limit more directly.To derive an evolution equation for the average concentration from the master equation would require evaluating the first moment of the master equation, and considering mean-field approximations to remove nonlinearities in the reaction term [40].One of the difficulties of this approach is to combine both stochastic reactions and stochastic jumps in the master equation.
We refer the interested reader to Refs [42,43], in which such an approach is undertaken.The continuum limit leading to Eq. ( 28) is thus expected to represent the average behaviour of the stochastic model under a mean-field approximation where there are no long-range correlations and only weak nonlinearities [40,46].In particular, it may not hold well if the network possesses far connections over distances of the same order as the whole domain.The limit ∆ξ → 0 means that distances between connected nodes remain small with respect to the whole domain of the network.In practice, most spatial networks predominantly connect nodes in close spatial proximity because there is usually a cost associated with building connections.The evolution of the concentration of molecules n (x , t ) along the x -axis is found from Eq. ( 28) by using Eq. ( 22) and noting that For ease of notation, we now omit space and time arguments, and implicitly assume that all spatial arguments can be either functions of ξ or x , as appropriate, via the one-to-one mapping x (ξ).We have so that Eq. ( 28) becomes Therefore, the concentration of molecules n (x , t ) obeys the reaction-diffusion-advection equation where are a diffusivity and advective velocity in x -space, respectively, such that the total flux of signalling molecules at location x and time t is The space dependence of D (x ) and v (x ) arises due to the network's varying node density via the metric g , and due to the node connectivity via D and v .Using Eqs ( 27), (25), and the correspondence between continuous and discrete jump probabilities, we can elucidate how the local diffusivity D (x i ) and advective velocity v (x i ) at node i depend on the connections of node i with other nodes and the local node density ρ(x i ) as Equations ( 31), (33), and (34) provide direct relationships between macro-scale reactiondiffusion-advection properties of signalling molecules along x and local network topography, which includes both placement in space of nodes of the network via the node density ρ(x ), and node connectivity via the first and second moments of jump probabilities.The absorbing boundary conditions in Eq. ( 7) assumed for the stochastic and deterministic discrete models carry over in the continuum limit, such that Initial conditions assumed in the stochastic and deterministic discrete models determine initial conditions η(ξ, 0) and n(x , 0) for the concentration of molecules via Eq.( 21).
The flux of molecules J ± (t ) reaching the boundaries of the domain at time t in the continuum limit has no contribution from the advective flux v n due to the boundary conditions (36).It is due to the diffusive part of the flux only, i.e., The sign of the flux J − (x , t ) at the left boundary of the continuum model is chosen to ensure consistency with the non-negative flux at the left boundary in the stochastic and determinisitc discrete models, see Eqs (8).
When all the jump probabilities q s i are set to zero except for s = 0, ±1, ±2, we obtain jumps to nearest neighbours and next-nearest neighbours only.In this case, using the fact that To ensure that advection is of the same order as diffusion in the continuum limit, the jump probabilities must be chosen such that the following biases are proportional to ∆ξ:

Numerical discretisation
To solve the reaction-diffusion-advection Eq. ( 31) numerically, we use a simple finite difference scheme based on a regular discretisation of x , M , and a regular discretisation of t , t γ = γδt , γ = 0, 1, 2, . ... The numerical solution n (x α , t γ ) is then stepped in time using a simple explicit forward Euler scheme with a first-order upwind approximation of the advective term , and a second-order central difference approximation of the diffusion term ∂ x [47].We ensure that δt and δx fulfill the Courant-Friedrichs-Lewy condition δt ≤ min δx

Results and discussion
We present a range of numerical results for the stochastic and deterministic discrete models and their continuum limit.We first validate the models on a regular network with space-independent jump probabilities prior to increasing the complexity of the considered numerical examples.In Sections 3.1 and 3.2, space and time units are considered to be arbitrary and omitted in the numerical values presented.In Section 3.3, our model is applied to signalling in the osteocyte network and we provide relevant space and time scales.

Space-independent jump probabilities
In this section we assume that every network node x i , i = 1, . . ., K is connected to its immediate neighbours only and that jump probabilities are space-independent.Therefore, q −1 i = q −1 , q 0 i = q 0 and q 1 i = q 1 for i = 2, . . ., K − 1, with all remaining q s i = 0.In addition, we assume that K is odd and the network is symmetrical about the central node ĩ = K +1 2 at x (K +1)/2 = 0.All molecules are initialised at the central node of the network at time t = 0: where N tot is the number of initial molecules and δ i ĩ is the Kronecker delta.These initial conditions allow us to consider the response of the network to an impulse signal generated at the centre of the network at time t = 0.In the continuum limit in which ∆ξ → 0, the initial condition of the reaction-diffusion-advection equation given by Eq. ( 41) becomes a Dirac delta function: For the finite difference scheme these initial conditions become N tot /δx at the discretisation point at x = 0, and zero elsewhere, where δx is the spatial discretisation and the number of discretisation intervals M is even.
The above assumptions result in constant advection v and constant diffusion D in ξ-space.In the case that there are no reactions (F = 0), Eq. ( 28) simplifies to a space-homogeneous diffusion-advection equation for η(ξ, t ).The response to an impulse signal of N tot molecules in the centre of an infinite domain in ξ-space is therefore given by the fundamental Gaussian solution This gives the following infinite-space exact solution n ∞ (x , t ) for Eq. ( 31) when nodes are positioned arbitrarily, the jump probabilities are space-independent, and there are no reactions: We first consider a regular network with We set jump probabilities as q 0 i = 1/3, q −1 i = 1/3−a ∆ξ, and q 1 i = 1/3+a ∆ξ for i = 2, . . ., K −1, with all other q s i = 0.The parameter a biases jumps to the right (or left) when positive (or negative), or jumps are symmetrical (unbiased) when a = 0.In this case the advection velocity v and diffusivity D are and we emphasise that they are of the same order in ∆ξ and ∆t .Figure 2 presents the concentration of molecules for the stochastic, deterministic, and continuum models for the two cases of unbiased and biased jumps.There is an excellent match between the deterministic discrete model and finite difference solution to the continuum equation, with variability of the stochastic model around these solutions.The infinite-space exact solution n ∞ matches the finite difference numerical solution until molecules begin reaching the boundaries and start being removed from the finite system.As expected, biased jumps result in molecules on average moving to the right.
Figure 3a illustrates the metric, node density, advection, and diffusion functions for this network.Figure 3b shows results for the evolution of concentration of signalling molecules within this network starting from the impulse initial condition in Eqs ( 41), (42).There is again good agreement between the models, with the exception of the infinite-space exact solution n ∞ that does not satisfy the zero Dirichlet boundary conditions on the finite domain.The time snapshots clearly demonstrate how the molecule concentration increases in regions with higher node density.This is expected (see Eqs ( 21) and ( 43)), and results in an increase in signal strength as the impulse propagates away from the centre of the network.This phenomenon and its relevance to the osteocyte network in bone is explored in more detail in Section 3.3 below.We now consider a network where the node density changes periodically.We choose where A > 0 and k > 0, and ξ(x ) can be calculated by inverting x (ξ) numerically.For the mapping x (ξ) to be bĳective, we require Ak < 1.As above, q 3 for i = 2, . . ., K −1. Figure 4a illustrates the functions g (x ) and ρ(x ) as well as the advection and diffusion functions determined numerically from the following expressions derived from Eqs (32): and Figure 4b presents a time snapshot of the molecule concentration of the stochastic, deterministic, finite-difference, and infinite-space solutions, and also presents the total number of molecules Concentration of molecules corresponding to the deterministic (gray), stochastic (blue), finite deference (FD, green), and infinite-space exact (red) solutions at t = 2, and conservation of molecules for these models over time.Parameters: 3 for internal nodes, for each model over time.As in Fig. 3, we see the influence of the node density on the signal concentration.Figure 4b shows that the finite difference method does not accurately solve the differential equation for this periodic network.This is due to the metric g (x ) being small away from the boundaries of the computational domain, which results in a nearly singular reactiondiffusion-advection equation and consequently a significant numerical loss of molecules in this non-conservative finite-difference discretisation scheme.A more sophisticated, conservative numerical scheme could address this issue.However, we emphasise that the deterministic discrete model is a fast and effective approach for solving the reaction-diffusion-advection equation that is not affected by small values of g (x ) and that remains exactly conservative.The small loss of molecules seen in Figure 4b in the stochastic and deterministic discrete models is due to molecules exiting the network at its left and right boundaries.

Space-dependent jump probabilities
In this section we consider space-dependent jump probabilities and determine how these probabilities can be selected to give rise to given diffusivity and advection in the continuum limit.This selection is not unique in general so for simplicity, we restrict to nearest neighbour jumps so that only q −1 i , q 1 i , and q 0 i = 1 − q −1 i − q 1 i are non-zero.In this case, we can rearrange Eqs (33), and ( 34) to obtain space-dependent jump probabilities q −1 i , q 1 i , and q 0 i as functions of D and v : where the right-hand sides are evaluated at x i .Appropriate bounds on ∆t and ∆ξ are required such that all the jump probabilities are non-negative.
For illustration, we solve the spatially homogeneous diffusion equation on an irregular grid using these expressions for jump probabilities.Assuming zero advection (v = 0) and constant diffusivity (D = D 0 ), Eqs ( 49)-(51) become For the periodic network given above in Eq. ( 46), the jump probabilities that give the diffusion equation with diffusivity D 0 in the continuum limit are , , We define the network over the interval [x 1 , x K ] = [−L , L ] and continue to use absorbing boundary conditions as in Eq. (36).We use the initial molecule concentration and round N i (0) = n (x i , 0)∆x i to the nearest whole molecule at each node x i for the stochastic model.This choice of initial condition enables comparison of the stochastic, deterministic and finite difference models with an exact solution on the finite domain.This exact solution to the diffusion equation with diffusivity Figure 5a presents a time snapshot of the deterministic, stochastic, finite difference and exact solutions for n (x , t ).As expected, the finite difference solution closely matches the exact solution.However, the deterministic and stochastic models deviate from the exact solutions where the node density ρ(x ) is high.In Fig. 5b we present the absolute error between the exact and deterministic solutions for K = 51 and K = 101 at time t = 2.With more nodes in the network, spacings between nodes become smaller and the discrepancy between the exact and deterministic solutions decreases from a maximum of 17% relative error with K = 51 to a maximum of 5% relative error with K = 101, consistently with the reaction-diffusion-advection equation being the continuum limit of the deterministic model.Figure 5c presents the variance Var n stoch i (t ) at time t = 2 calculated using 100 stochastic realisations of this model, each with N tot = 1000 molecules.This is shown alongside the expected variance Var (n (t )) that is computed with the right-hand expression of Eq. ( 13) using the exact solution n (x , t ) from Eq. ( 55) in place of n disc i (t ).It is clear from Fig. 5c that the variance increases in areas of the network with high node density.This result may not be intuitive; it emphasises the importance of considering stochastic models, particularly in biological contexts where molecule numbers are small and fluctuations can be large.Mathematically, this result is due to the uneven distribution of nodes in space and the fact that the variance of number of molecules at a node is quadratic in P i (t ) in Eq. (11).The deterministic and continuum solutions cannot capture this node-density-dependent behaviour of molecular fluctuations.
Yates et al. [49] showed using first passage time (FPT) problems that in order to achieve constant diffusivity in the continuum limit on an irregular network, left and right jump probabilities should be selected as where To show that the first passage time jump probabilities q ±1 i in Eqs (56) converge to q ±1 i in Eqs ( 52) in the limit ∆ξ → 0, we first note that given h i = x i − x i −1 , we have Substituting Eqs ( 57)-( 59) into Eqs ( 56), expanding g , we obtain, in the limit ∆ξ → 0, which matches the expressions q ± i in Eqs (52).We have computed the FPT solution, n FPT , using the jump probabilities from Yates et al. [49] corresponding to our periodic network.Figure 5d presents the disparity between the exact solution of Eq. ( 55) and the FPT solution for K = 51 and K = 101 at time t = 2.As for our deterministic model, the error decreases with more nodes in the network.The error in our deterministic model is significantly larger than the error using the jump probabilities of Yates et al. [49].However, we emphasise that our approach provides explicit expressions for q −1 i , q 0 i and i in terms of an arbitrary desired advective velocity v (x ) and diffusivity D (x ) on an arbitrary network (Eqs ( 51)-( 50)).

Signal propagation mechanisms in the osteocyte network
Signals generated by osteocytes within bone tissue propagate through the osteocyte network to the bone surface, where they induce bone formation or bone resorption [32,33,34].In this section we propose a simple, illustrative model of signal generation and propagation through a one-dimensional osteocyte network (Fig. 6).We investigate how signals arriving at the boundaries of the network are influenced by local perturbations occurring at internal nodes (osteocytes) of the network.The one-dimensional network that we consider may represent osteocytes distributed across the cortical wall of a long bone such as a mouse tibia, between the bone surface of the marrow cavity (endosteal surface) and the outer bone surface (periosteal surface).In the mouse tibia, osteocyte density is greater near the endosteal and periosteal bone surfaces and lowest in the middle of the cortical wall, which is approximately 200 µm thick [35,36].
We choose the x -axis such that the origin represents the middle of the cortical wall, and consider where L = 50 µm is a spatial scale parameter.With this choice the density of osteocytes increases away from the origin as ρ(x ) = e |x |/L /(L ∆ξ).This choice is not intended to be an accurate representation of osteocyte density measurements in mouse tibia [36], but it allows us to qualitatively capture density heterogeneities in real bone while allowing closed-form expressions for g (x ) and therefore D (x ) and v (x ), as well as the long-time behaviour of the solution in infinite space.Osteocytes are highly connected to neighbouring osteocytes [10].We set q 0 i = 0.36, q ±1 i = 0.16, q ±2 i = 0.16 for all internal nodes, with all remaining q s i = 0.This means that each osteocyte is connected to its nearest and next-nearest neighbours, and there is no jump-probability bias.
The metric, node density, advection and diffusion functions for this network are presented in Fig. 7.The jump probabilities and discretisation parameters of the stochastic and deterministic models are chosen so that advective velocities v (x ) are of the same order of magnitude as those arising in the fluid flow model of van Tol et al. (about 20 µm/s) [35], and diffusive root mean square displacements 2D (x )t over one second are of the order of 20-80 µm.While these values are commensurate with effective fluid flow transport properties [35], other parameter choices could be made to represent other types of signalling through the osteocyte network such as calcium signals [50,51].We note here that the results we present in this section illustrate qualitative behaviour of the model that do not depend on the space and time scales chosen.By rescaling time, space, and densities in Eq. ( 31), it is possible to modify the scales of the axes in all the figures presented below [52].
The network's metric g (x ) and diffusivity D (x ) have a cusp at x = 0 and the advective velocity v (x ) is discontinuous at x = 0.Because of these singularities in D (x ) and v (x ), the finite difference solutions of the reaction-diffusion-advection equation are inaccurate compared to those provided by the deterministic model (see Section 3.1).We therefore omit numerical results from the finite difference simulations in this section.
We first analyse the robustness of signal propagation under perturbations of the network.We assume no chemical reactions (F = 0) and consider that an impulse signal is generated at time t = 0 at the centre of the network (x = 0).Figure 8 shows schematics of the considered network perturbations (Fig. 8a) and corresponding snapshots of the molecule concentration at time t = 2.5 (Fig. 8b). Figure 8c shows the time signature of the fluxes J − (t ) and J + (t ) of molecules reaching the left and right boundary nodes in response to the initial impulse.These responses represent signals that could induce bone formation or bone resorption at the bone surface.The left column of Fig. 8 relates to the unperturbed network.The initial impulse splits into two peaks propagating left and right toward the network boundaries (Fig. 8b).While diffusion tends to disperse molecule concentration peaks (cf.Fig. 2), here the increasing node density encountered by the molecules as they propagate away from the origin counters their dispersion, and creates well-defined outward-moving peaks evolving into a travelling wave of constant shape and height on an infinite domain (see Appendix A).The centre column of Fig. 8 relates to a node i = 20 being removed in the left half of the network.Signalling molecules are still transmitted through to the left boundary i = 1 due the presence of next-nearest neighbour connections between nodes i = 21 and i = 19.However, molecules tend to pile up at node i = 21 because they have fewer opportunities to leave this node toward the left, resulting in an increased concentration right of the node i = 21 and a decreased concentration left of the node i = 19.This results in a lower peak in the flux signature J − (t ) at the left boundary (Fig. 8c).The right column of Fig. 8 relates to second-neighbour connections of node i = 20 being severed.In this case molecules progressing toward the left boundary pile up before node i = 20 since they transfer over the severed region less easily, resulting in an even lower flux J − (t ) at the left boundary.These local changes in the network may represent micro-damage in bone tissue such as osteocyte death or micro-cracks severing connections.Micro-damage in bone is removed by a bone remodelling process initiated at the bone surface [53,33].Differences seen in our simulations between the signal received at the left and at the right boundaries due to an off-centre damaged osteocyte could thus enable localisation of the damage and help guide repair processes, such as triggering a bone remodelling event to be initiated closest to the damage.
We now investigate the influence of localised changes in molecule generation by the osteocyte network, and how these changes affect the signals received at the bone surface (network boundaries).We assume a steady, homeostatic state, in which each osteocyte generates a background of signalling molecules with rate λ 1 (number per unit time).The molecules are also assumed to degrade with a first-order reaction rate λ 2 (probability per unit time), so that This homeostatic state models osteocytes subjected to normal mechanical stimulation.The constant creation and elimination of molecules as they diffuse through the network generates a fluctuating quasi-steady state in the stochastic model, with an average concentration of molecules along x that is affected by the node density and the absorbing boundary conditions, as well as the strength of reaction rates λ 1 and λ 2 (Fig. 9a).In the stochastic model, the fluxes J ± (t ) of signalling molecules at the network boundaries are highly fluctuating quantities in the presence of stochastic reactions because we assume only one reaction per node per time step ∆t .We therefore report the flux moving average values over a time window of 10 s.
If an osteocyte at the centre of the network produces an additional impulse in the number of molecules at a certain time, this impulse will split into two peaks propagating left and right toward the boundaries and generate a peak in J ± (t ) similarly to Fig. 8.However, the intensity and duration of this impulse response depends on the reaction rates λ 1 and λ 2 (Fig. 9).Indeed, the continual creation and degradation of molecules at each node acts as a chemical reservoir that tends to equilibrate the number of molecules at each node [40,54].Higher reaction rates lead to a faster return to equilibrium.This results in a more attenuated response (Fig. 9b), but also in less broadening of the signal in time (Fig. 9c).
Figure 9b shows that the intensity of the impulse response J + (t ) at the network boundary may be of similar order of magnitude as inherent fluctuations of J + (t ), and may therefore not be easily detected.In Fig. 10 we show several stochastic realisations of the signal response J + (t ) received at the right boundary of the network following either an impulse increase in molecules at the central node at time t = 50 s, or a step increase in the rate of molecule production λ 1 at the μ μ Figure 9 -Influence of reactions on signalling within a network with exponentially increasing node density for the three cases of no reactions (λ 1 = 0 and λ 2 = 0, blue curves), reaction rates of λ 1 = 2 s −1 and λ 2 = 0.005 s −1 (red curves), and reaction rates of λ 1 = 20 s −1 and λ 2 = 0.05 s −1 (green curves): (a) Concentration of molecules at steady state in the deterministic model; (b) Flux of the deterministic and stochastic models at the right boundary starting from the steady state and with an impulse of 600 molecules added at x = 0 at time t = 50 s.The curves are shifted by J + (0) so they start from the same baseline.(c) Flux of the deterministic model at the right boundary as in Fig. 9b, where the curves are shifted by J + (0) and scaled by J m + = max t J + (t ) − J + (0).(Parameters as given in Fig. 7. Stochastic fluxes are averaged over a time window of 10 s).central node for t ≥ 50 s, starting from a (fluctuating) initial steady state at t = 0.The first signal production represents an osteocyte generating a sudden, short-lived burst of signalling molecules, such as to apoptosis [34].The second signal production respresents a continued response in time, such as due to the detection of local micro-damage, or mechanical overstimulation.In the deterministic model, both types of signals result in a noticeable (average) response, but in the stochastic model, some stochastic fluctuations in steady state may be misinterpreted as a response to an impulse signal.In contrast, the response to the step increase in production rate generates a stronger signal-to-noise ratio in the stochastic model.

Conclusion
In this paper we have presented stochastic and deterministic versions of random walk models for the reaction and propagation of signalling molecules on one-dimensional spatial networks with arbitrary node placement and connectivity.We have derived a reaction-diffusion-advection differential equation that is the continuum limit for these models.Solutions of the reactiondiffusion-advection equation capture the average behaviour of the stochastic model well when spacings and connections between nodes are small, or equivalently, when the network is dense.
The continuum limit elucidates how signalling molecules propagating through spatial networks are transported in physical space through diffusion and advection properties that depend on local node density and connectivity within the network.The relationships we derived between spatial network properties and diffusion and advection can help us understand average transport properties in physical communication networks.Our results highlight that the concentration of molecules in space and their diffusivity depend strongly on the node density of the network.Node density can thus have important consequences for the development of spatial patterns of morphogens in addition to the influence of network connectivity [29].Node density impacts the efficiency with which molecular factors diffuse and establish in biological tissues to perform their function.As already mentioned in the introduction, in plants, for example, the functional behaviour of the network of cell-to-cell communication controls average tissue-scale behaviours such as tissue growth and defense mechanisms [16], as well as patterning and organ development [7].In distributed such as cardiovascular, respiratory and root systems, spatial network properties will influence average propagation times of nutrients and waste products [14].
In the example of bone tissue, the osteocyte network regulates the mineralisation of bone, which provides compressive strength and fast calcium storage and retrieval [55,56].Experimental data on osteocyte networks and mineral density exhibit patterns suggesting complex diffusion and reaction behaviours of mineral precursors and inhibitors through the osteocyte network [57], that could lead to future analyses based on our modelling approach.
While the deterministic discrete model and continuum reaction-diffusion-advection equation provide accurate representations of the average behaviour of the stochastic random walk model, they do not account for fluctuations in molecule concentrations inherent to stochasticity in both transport and chemical reactions.The intensity of fluctuations was found to depend strongly on the network's node density.This could affect signal-to-noise ratios significantly in signals propagating through cellular networks, in which a baseline of stochasticity can be large.Some chemical reaction behaviours such as regulatory switches and Turing patterns can be triggered by fluctuations [58,59].The consideration of the spatial arrangement of nodes of the network in stochastic signal propagation in these situations could lead to new behaviours that are of high interest.
Numerically solving inhomogeneous reaction-diffusion-advection equations with rapidly varying diffusivity and advective velocity can be challenging.Sophisticated conservative finitevolume methods or finite element methods are often used for this purpose, since standard finite difference methods can fail to provide accurate solutions, as illustrated by our results when the metric is close to zero and the diffusivity and advective velocity have singularities.However, our results show that the deterministic discrete model can be used as an effective, conservative, and explicit numerical method for solving the reaction-diffusion-advection equation in all situations.Recent work has utilised a similar approach to solve nonlinear advection-diffusion equations on regular grids [41].Our approach enables us to obtain expressions for nearest-neighbour jump probabilities that will achieve desired inhomogenous advective velocity and diffusivity in the continuum equation for a given, arbitrary grid.
We illustrated a particular application of our model to the osteocyte network in bone.In this system, positional information about where signals originate from deep inside the network in response to damage is important so that this damage can be repaired [53].Signals generated by osteocytes within bone tissue in response to damage propagate to the bone surface, where they induce repair through bone resorption and bone formation.The spatial location of these signals is critical, so that bone resorption is not initiated where bone tissue is needed most, and so that micro-damage within bone tissue induces targeted remodelling.Our investigation of networks in which internal nodes are perturbed shows that signals received at the boundaries of the network can indicate the location of the perturbation.However, these signals may also be difficult to detect depending on the nature of signal generation and the presence of chemical reactions.The considered reactions contribute significantly to noise in the signal received at the network boundaries, with noise increasing with reaction rate.The noise can result in an impulse of signalling molecules at the centre of the network being difficult to detect at the network boundaries.However, a step signal, for example an increase in the reaction rate at the central node of the network, generates a stronger signal-to-noise ratio at the network boundaries.Our model of the osteocyte network is currently limited by only accounting for one spatial dimension.In a one-dimensional network, localised damage is always necessarily traversed by propagating signals.In higher dimensions, localised damage may affect propagating signals less as there are more paths available to bypass the damage.Experimental datasets of 3D osteocyte network architectures from confocal microscopy have been employed to simulate fluid flow through the lacuno-canalicular network in which osteocytes reside [35,36].While our model considers a simplified, one-dimensional situation, it still represents important mechanisms by which the osteocyte network is believed to operate in real bone tissues.Some of these mechanisms have not been previously considered, such as the influence of disruptions to the network, chemical reactions, and stochasticity.
The extension of our model to higher dimensions and to dynamic, adaptive networks [18] is of strong interest as it would enable more comprehensive modelling of several network systems.Deriving continuum limits for arbitrary networks in higher dimensions is a nontrivial extension but it can be anticipated that average transport properties of signals may also be governed by reaction-diffusion-advection equations, since these types of equations express conservation laws that hold at the discrete level [60,61].In bone, such models may be able to further our understanding of how signals propagating through the osteocyte network contribute to bone formation and resorption processes, and how osteocytes set the mechanical memory of bone [62].
Further avenues for future work include the consideration of more detailed transport problems, such as propagation of signalling molecules along the connections between the network's nodes, which result in partial differential equations on networks [18].

Data accessibility
Key algorithms used to generate results are freely available on Github at https://github.com/prbuen/Mehrpooya2023_RandomWalk1DNetwork.The wave is generated after an impulse initial condition at the origin, and it travels toward the right with speed v (t ) = L /(2t ).Nodes of the network are distributed according to Eqs (60).
Thus, the exponential increase of node density with x in this network counters precisely the diffusion-driven dispersion of the initial impulse, and gives rise to a travelling wave of constant profile, but decreasing speed.

Figure 1 -
Figure 1 -(a) Schematic of a random walk on irregularly spaced nodes x i , with only the incoming connections to node i displayed.(b) The regular network of nodes ξ i is mapped to the x i network by x (ξ).

Figure 4 -
Figure 4 -Network with periodic node density, x (ξ) = ξ + A sin(k ξ): (a) Metric, density, advection and diffusion functions.Network node positions are marked on the x -axis.(b)Concentration of molecules corresponding to the deterministic (gray), stochastic (blue), finite deference (FD, green), and infinite-space exact (red) solutions at t = 2, and conservation of molecules for these models over time.Parameters: q 0

Figure 6 -
Figure 6 -Schematic osteocyte network across the cortical wall of a mouse tibia cross-section.(a) Cortical wall of a mouse tibia cross section (grey) showing the x -axis, signalling molecule fluxes J − (t ) and J + (t ) at the endosteal and periosteal bone surfaces, respectively, and the location of osteocyte network damage (red cross).(b) Osteocytes along x connected into a one-dimensional network.The osteocyte in blue represents a stimulated osteocyte generating signalling molecules transfering to neighbouring osteocytes by the jump process.The osteocyte and connection in red represent localised network damage.

Figure 8 -
Figure 8 -Perturbations of a network with exponentially increasing node density: (a) Schematic of the network depicting the unperturbed case, removal of node i = 20, and removal of secondneighbour connections at i = 20.(b) Concentration of molecules corresponding to the deterministic (gray) and stochastic (blue) solutions at time t = 2.5 s for the three networks from an initial impulse with N tot = 10, 000 molecules at x = 0.The presented infinite-space exact solution (red) in all three columns relates to the unperturbed network.(c) Flux corresponding to the deterministic solution at left boundary (gray line) and at right boundary (black dashed line) for the three networks.(Parameters as given in Fig. 7).

Figure 10 -
Figure 10 -Response to impulse and step signals from steady state within a network with exponentially increasing node density: flux of the deterministic and stochastic models at the right boundary for an impulse of 600 molecules added at x = 0 at t = 50 s (dark gray and light gray), and a step increase in λ 1 from 2 s −1 to 20 s −1 at x = 0 for t ≥ 50 s (dark green and cyan).Each case includes two stochastic realisations.The dashed black line indicates the mean of the stochastic flux at the steady state, and the blue dashed lines indicate one standard deviation from the mean.(Parameters as given in Fig.7with λ 1 = 2 s −1 and λ 2 = 0.005 s −1 for all nodes except at x = 0, and a time averaging window for the stochastic flux of 10 s.)

Figure 11 -
Figure 11 -Long-time wave profile φ(y ) given by Eq. (A 64) with N tot = 10, 000 and L = 50 µm.The wave is generated after an impulse initial condition at the origin, and it travels toward the right with speed v (t ) = L /(2t ).Nodes of the network are distributed according to Eqs(60).