Biological network growth in complex environments: A computational framework.

Spatial biological networks are abundant on all scales of life, from single cells to ecosystems, and perform various important functions including signal transmission and nutrient transport. These biological functions depend on the architecture of the network, which emerges as the result of a dynamic, feedback-driven developmental process. While cell behavior during growth can be genetically encoded, the resulting network structure depends on spatial constraints and tissue architecture. Since network growth is often difficult to observe experimentally, computer simulations can help to understand how local cell behavior determines the resulting network architecture. We present here a computational framework based on directional statistics to model network formation in space and time under arbitrary spatial constraints. Growth is described as a biased correlated random walk where direction and branching depend on the local environmental conditions and constraints, which are presented as 3D multilayer grid. To demonstrate the application of our tool, we perform growth simulations of a dense network between cells and compare the results to experimental data from osteocyte networks in bone. Our generic framework might help to better understand how network patterns depend on spatial constraints, or to identify the biological cause of deviations from healthy network function.

Introduction

Complex biological networks such as the neural connectome are striking examples of large-scale functional structures arising from a locally controlled growth process [1, 2]. The resulting network architecture is not only genetically determined, but depends on biological and physical interactions with the microenvironment during the growth process [3, 4]. In this context, evolution has shaped diverse spatial networks on all length scales, from the cytoskeletal network in cells [5], to multicellular networks such as the vascular system [6] or the osteocyte lacuno-canalicular network [7], to macroscopic networks of slime molds [8], mycelia [9] and plants [10]. Sophisticated imaging techniques together with large-scale automated analysis provide increasingly detailed views of the architecture of biological networks, revealing e.g. how neurons are wired together in the brain [11]. After extracting the topological connectivity from such image data, quantitative methods from the physics of complex networks can be applied to compare different types of networks and to uncover common organizational principles [12-15]. To understand the functional role of networks in the context of evolution, however, it is not sufficient to characterize static network structure [1]: we need to be able to infer the dynamics of the underlying biological growth process that defines this structure.

This poses a major challenge, since the local dynamics of the growth process is in most cases not experimentally accessible, except for simple model systems. Here, computer simulations provide a solution: they allow to connect the observed patterns to the underlying process by performing in-silico experiments. Different hypotheses can be tested by systematically varying the local growth rules in the simulation and analyzing the resulting network patterns. Such a computational approach also helps to overcome an important limitation of network physics with respect to spatial networks: the properties of complex networks are often calculated with respect to canonical random graphs such as the Erdös-Renyi [16] or Watts-Strogatz model [17]. These models are defined by the topological structure in an adjacency matrix, but usually neglect spatial constraints [18, 19]. A generic model of network growth under spatial constraints would thus address two important problems: enabling a meaningful quantification of spatial network patterns, and linking the observed patterns to the underlying biological growth process.

The first computer simulation of spatial network development, published in 1967 [20], already included the role of active and repulsive cues. Since then, several more advanced approaches to model spatial network growth have been developed. One type of model implements network growth by formation of edges between pre-existing nodes in space [21, 22] and was e.g. used to develop large-scale models of branching neurons [23] or leaf vascular networks [24]. Another type of approach is based on growth processes emanating from predefined seed points with branching and merging and has been used to model neural network development [25, 26], fibrous materials [27] or the development of branching organ structures [28]. Most of these existing models are targeted towards a specific domain such as neural networks and synapse formation, or include only simplified interactions with the environment. The current state of modeling spatial network growth was recently summarized in [1], arguing that there is “a crucial lack of theoretical models”.

The aim of this work is to develop a generic framework for biological network development in space and time under arbitrary spatial constraints. We propose a probabilistic agent-based model to describe individual growth processes as biased, correlated random motion with rules for branching, bifurcation, merging and termination. We apply mathematical concepts from directional statistics to describe the influence of external cues on the direction of individual growth processes without restricting the model by making too specific assumptions about these cues. Structural and geometric factors are obtained directly from real image data rather than formulated explicitly. The model is implemented as a computational framework that allows us to perform simulations for diverse types of biological networks on all scales, and to monitor the evolving spatial structure and connectivity.

Results

Growth process

The elementary process of biological network formation is the outgrowth of new edges in space from existing nodes, such as cells. Each individual growth agent, for example a neuronal growth cone, randomly explores space and senses soluble signals as well as physical cues and obstacles. This process can be formulated as a biased correlated random walk [29]: the movement of edges is biased by attractive or repulsive signals and structural cues, and correlated due to directional persistence of edge growth. Edges turn into trees if branching and/or bifurcations occur, such as in growing neurons, osteocytes or sprouting angiogenesis. Finally, if edges can join other edges and form new nodes (e.g. synapses), a connected network emerges. The elements of this discretized, agent-based formalization of spatial network development are summarized in Fig 1. In the following, we describe our mathematical approach for modeling the growth process.

Fig 1

Discrete local events during spatial network development.

Two growth agents (green and blue, e.g. growth cones, filopodia, cells) start at time t and randomly explore their environment in a biased correlated random walk. The environment is represented as a discrete grid (grey lines), while the network grows in continuous space. Blue and green dots correspond to subsequent positions of the growth cones. Cones can branch off new daughter cones, bifurcate/split in two new cones, merge with existing edges, or terminate. These events leads to formation of new nodes in the network, shown as red dots. Over time, the individual processes form a connected network of nodes and edges in space.

In general, such a stochastic motion under the influence of external forces is described by a partial differential equation (PDE), the Fokker-Planck equation [30] corresponding to the Langevin equation [31]:

Eq 1 is a drift-diffusion equation of motion where the k-dimensional vector is a drift, and D the anisotropic and diagonal, k × k- dimensional diffusion tensor. This equation describes the time development of a probability density function of diffusing particles under the influence of a drift and diffusion . If drift and diffusion are both constant in time and space, and if the initial condition is the k-dimensional Dirac delta function, then Eq 1 solves to

Eq 2 describes the time development of a random movement for constant drift and isotropic diffusion starting at initial position . In a complex environment, however, both drift and diffusion can depend on time and position, and there is no simple analytic solution for the time development. In this case, a numerical approach can be taken by sampling individual movement traces as non-Markovian correlated processes [32]. The drift-diffusion movement is then discretized into a stepwise process where for each growth step, a new direction is drawn from a probability density function (PDF). For constant drift and diffusion (Eq 2), the corresponding PDF, or growth kernel, is a Multivariate Gaussian distribution (MVG) with mean vector , covariance Σ = D t and determinant . In the general case, assuming that the spatio-temporal change for and is small between two consecutive time steps, every discrete growth step has a different time- and position dependent PDF. The fundamental solutions Eq 2 of the PDE (Eq 1), or Eq 3 for non-diagonal covariance, can thus be used as random kernels in each step from which the subsequent growth directions are drawn. Their and Σ can be expanded so that they encode for all spatial signals and constraints. A more detailed derivation of this approach starting from Eq 1 is provided as supporting information in S1 Text.

Substrate and signaling cues

The individual PDFs that describe all the different internal and external constraints and cues that act on an edge at position are formulated using the Multivariate Gaussian Distribution (MGD) as kernel:

The MGD is a probability density function from which random vectors of dimension k are drawn. Such a vector is shifted by the mean vector and anisotropically distributed. The anisotropic distribution of directions is described by the metric of the k × k covariance matrix Σ with determinant |Σ| and for its eigendecomposition as real symmetric matrix Σ = A− ΣA with eigenvectors and eigenvalues Σ = diag(λ1 < ⋯ < λ).

Depending on the edge’s current position , the mean vector points to the most likely direction the edge will be moved in the subsequent step, while the covariance Σ reflects the impact of the local structure onto its movement in the complex environment.

Following the mathematical rules for distributions, the PDF of a random vector that is the sum of independent random vectors is the convolution of the PDFs of the individual random vectors:

All local, external signals and structural properties that are described by individual PDFs in the form of a MGD can thus be combined into a single PDF that still has the form of a MGD [33]. As a result, a new growth direction that includes the contribution of all cues and constrains can be drawn from this combined MGD with mean vector and covariance

These parameter pairs can represent different external cues, the geometric constraints of the substrate, as well as edge-specific properties. In Table 1, we give some examples for external and internal cues and the structure of the environment that can be represented by such “Mean Vector—Covariance” pairs.

Table 1

Mean vector—Covariance pairs.

External Cue	μ→	Σ
Drift	D→(p→(x,y,z))	σx,y,z → 0
Signaling molecules	±∑mn(p→m,is-p→m,i-1s)n	∑mnΣ_m
Internal Cue	μ→	Σ
RW - Random Walk	0→	diagonal, σx,y,z > 0
P - Persistence	∑i=mnp→_i	isotropic, diagonal, σx,y,z > 0
Structural Cue	μ→	Σ
Structure tensors	∇→p(x,y,z)	Σ=αR^TΣdiagdualRΣ_diag=RΣ_i,jR^TΣ_i,j=δ_ipδ_jpi,j∈(x,y,z)

Table 1: Examples how external, internal and structural cues interact with a growing edge at position in 3D. The physical properties are interpreted as multivariate Gaussian distributions with a mean vector and covariance Σ. All signals are combined into a single PDF by convolution. We distinguish between external signals such as drift and signaling cues, internal properties of an edge such as its random motion and stiffness, and the structural properties of the surroundings. A more detailed version is provided in S2 Table.

Typical examples of external cues are drift vector fields and ensembles of signaling molecules. For the case of a constant external vector field, the mean vector points towards the local field orientation, while its length scales with the field’s strength. Its corresponding covariance is symmetric, proportional to the unit matrix , and the entries tend to zero (σ1,…, → 0). The MGD is then a k-dimensional delta distribution that biases the growing edge in the direction of the vector field.

Internal signals describe the properties of the growing edge, reflecting its characteristic random motion and its stiffness or persistence. The physical interaction with the substrate is encoded in the structure tensor and the local derivatives (gradients) of the surroundings. The structure tensor takes the place of the diffusive tensor in the MGD, while the mean vector is the local gradient . This concept, which is widely used in image analysis (e.g. anisotropic diffusion) [34-36], allows us to include arbitrary structural information by representing it as an image volume and calculating the corresponding structure tensor.

The combination of the different PDFs through convolution (Eq 4) results in a single MGD from which the subsequent growth direction is drawn. This combined PDF distributes random vectors around a mean vector with a length that is the sum of all mean vectors of the individual signals. The shape of this shifted multivariate normal is expressed by its covariance.

Step length

The Multivariate Gaussian Distribution in Eq 3 covers the full , but sampling the entire space to draw the next growth direction would not only be prohibitively slow, but would also generate unphysical step size distributions. This problem can be solved by integrating over the radius to generate a sampling of the PDF from a subset of (the unit sphere) while retaining the full spatial information contained in its covariance and mean. An additional advantage is that growth direction and step length are separated into two probability distributions which can then be treated independently.

Distributions on the unit sphere can be obtained by embedding in . For a random vector , this could be its conditional distribution on , or its angular projection . For the family of Fisher-Bingham distributions [37], the Bingham distribution [38] is the conditional distribution, while the Kent distribution is the projection of the general Fisher-Bingham distribution [39]. Although these distributions are well studied and are commonly used for data fitting and directional statistics [40], the Kent distribution is not suitable for simulations, as it requires a rejection method [41]. Instead, the MGD can be transformed to spherical coordinates, which allows us to calculate a spherical marginal distribution MGDangular by integrating over the radial component to obtain a projection of the MGD. Starting from Eq 3, the k-dimensional coordinate basis is transformed to with , and . The transformation and integration of Eq 3 results in the following angular MGD:

Carrying out the integration yields the general angular Gaussian distribution (GAGD), introduced and explicitly calculated for dimensions k = 1, 2, 3 in [42], which is the angular marginal distribution to the MGD that projects the directional probability information onto the -sphere: The last term of Eq 8 is a function For the 3-dimensional case k = 3, , this results in where ϕ(.) is the standard normal PDF and Φ(.) the corresponding CDF, the cumulative distribution function. For additional detail, please see S1 Text.

This definition separates the radial distribution from the angular distribution. While we use a constant step length in our simulations, this separation would also allow for variable step length PDFs without losing any information of the spatial distribution of external cues with large mean vectors . In Fig 2, we show three distributions on the unit sphere for three dimensions with different mean vector and covariance Σ as examples for shifted isotropic signals (persistence), highly anisotropic covariances due to spatial constraints such as edges and walls, and peaked distributions for guidance cues towards specific directions. The structural information is projected onto the unit sphere while retaining the full information about Σ and μ, as shown in Eq 8—e.g., a bimodal MGD still retains a bimodal marginal MGDangular due to its covariance.

Fig 2

Probabilistic framework to determine the next growth direction.

Probability distributions are modeled as multivariate Gaussian distributions (MGD), shown in a), with mean and covariance determined from the discrete simulation grid that describes structural and soluble cues. Individual MGDs are combined by convolution, transformed to spherical coordinates, and projected onto the unit sphere (bottom). Examples shown on the right hand side include structural guidance along a plane (e.g. a tissue boundary, b), a unidirectional drift towards the viewer (e.g. a growth factor gradient, c), and persistence due to memory of past growth directions (e.g. bending stiffness, d). The three probability distributions are merged by convolution, restricted by the aperture of the growth cone, and sampled on the corresponding segment of the sphere (e). The next growth direction is drawn from this combined, restricted distribution.

After motivating the mathematical kernel for the growth process and the interaction with the substrate, we now introduce the rules that enable the growing edges to bifurcate, branch off daughter processes, to form new connections (e.g. synapses) with other edges, and to be persistent and locally self-avoiding.

Persistence

During edge outgrowth, the finite bending stiffness of edges causes the growth direction to be correlated to the previous directions, resulting in persistence of the movement. In absence of a drift vector , and since the eigenvectors of the covariance matrix Σ are independent, the random growth process is a Wiener and Markov Process, and still remains a general interpretation of a Wiener process for . Persistence of the growth direction is equivalent to introducing a drift that is the sum of the last t growth directions . This results in an effective bending stiffness of the growing edges, and its persistence length is expressed by the Kuhn length [43], where δl is the length between two steps and δl〈cosθ〉 is the projected length at i towards j (Fig 1). This directional correlation does however not imply self-avoidance within the persistence length scale, as there still is a non-zero probability for an edge to loop back and collide with itself within the persistence length.

Self-avoidance

Local self-avoidance is the phenomenon that growing edges in a biological network, e.g. outgrowing neurites or branching vessels, do not loop back onto themselves [44]. In our growth model, local self-avoidance is achieved by restricting possible growth directions to a solid angle Ω around . The aperture described by Ω is equivalent to the maximum allowed curvature of a growing edge. This means that the heuristic growth process is no longer a Markov process [45], as it loses its reversibility, but now becomes self avoiding within the persistence length. Fig 2 shows the effect of persistence on the angular distribution MGDangular and how the aperture is restricted by Ω.

Connectivity

New nodes in a growing network form by branching or bifurcation of growing edges, or when a growing edge merges with another edge. We distinguish between “branching”, the branching-off of a new edge from a growing edge where the parent edge maintains its growth direction, and “bifurcation”, the splitting of a growing edge into two daughter edges growing in two different directions. Both can occur at every growth step with a probability that depends on edge properties (e.g. age) as well as external cues. Whenever an edge is within a certain distance to another edge, it can merge and form a new node, again with a probability that can depend on edge properties and external signals. Finally, growing edges can terminate and form an end-node, again with a probability that depends on external signals, such as obstacles or the density of already existing edges. Fig 1 schematically shows the time development of a growing network where edges branch several times, merge, and terminate.

Implementation

We next describe the computational implementation of the theoretical framework for spatial network growth described in the previous sections. Key input parameters for the simulation are the spatial environment and the seed points (e.g. cells) that define from where the first edges start to grow. The spatial environment is provided in the form of 3D image volumes that contain information about the physical structure of the growth environment (e.g. tissue structure), as well as about signals that influence the growth process (e.g. growth factors). This grid-based representation of the environment in an otherwise continuous model, as illustrated in Fig 1, makes it possible to study network growth in real tissue environments as measured e.g. with confocal or light sheet microscopy without explicitly stating their geometry by including them as scalar data on a 3D grid or “image” layer. Besides these image data and the seed locations and directions, the simulation requires a number of model parameters such as step length, aperture and degree of persistence as input, listed in Table 2. The aperture is the opening angle of the “field of view” of the growth cone. Smaller values limit the possible directions it can turn to, resulting in self-avoidance as described e.g. for neurons [44]. The angles for branching and bifurcations depend on the underlying biological mechanism, with large angles e.g. for vessel sprouting [46] or small angles for actin-crosslinker mediated branching of cell processes [47]. The corresponding probabilities for branching and bifurcation describe how likely these processes occur in a growing edge, which in biological systems can be both intrinsically determined (e.g. typical rates in osteocyte processes [14]) or externally controlled as in neurite outgrowth [2]. The parameter “memory” sets the number of previous steps that contribute to the next growth direction and determines the persistence of the direction. Biologically, persistent motion can arise either from physical constraints (e.g. persistence length of microtubules) or from the intrinsic persistence of the biological process, e.g. in cell migration. All parameters are set at the beginning of the simulation but can change for each growth cone as a function of time or external signals to capture the corresponding biological mechanism in the simulation.

Table 2

Simulation parameters.

Parameter	Typical Value	Biological Meaning
Aperture	Ω = 25 deg	Growth cone mobility, self avoidance
Step size	l∈R, l ∼ f(μ, σ)	Growth cone velocity
Splitting angle	θbran = 75 deg, θbif = 75 deg	Angle at which cones bifurcate or new branches form
Search field	3x3x3 cube	Range of interactions with the environment
Probabilities	pbran, pbif = 2%, pterm, preac = 0.1%	Probability to branch/bifurcate/terminate/reactivate
Memory	m = 5	Directional persistence of growth cone movement
Seed parameters	X, Y, Z, ϕ, θ	Seed positions and directions of growth processes

Table 2: List of parameters of the implementation with typical values, and their meaning in a biological context.

The PDFs for the growth processes are derived from the image volumes by calculating gradients and the structure tensor (Fig 3). These image volumes can contain multiple channels to describe not only structure, but also scalar and vector fields to describe forces or signaling particles. The spatial growth environment can also be dynamic, either via external time-dependent changes or via changes introduced by the growing network itself. Notably, it is by this reciprocal feedback between the growing network and its spatial environment that the network architecture turns into an emergent property that cannot be encoded in the growth rules alone.

Fig 3

Mapping the complex growth environment onto a multilayer simulation grid.

Structural and soluble cues and other signals are provided as 3D image volumes of real or simulated data (left). The six components of the structure tensor define the guidance cues from the substrate structure, whereas the three gradients of the images of soluble cues, e.g. growth factor concentrations, define the resulting directional bias. The feature sizes σ used for the filters determine the resolution of the cues. Other signals can be integrated by computing the corresponding features (filters) that define the resulting growth cue. The image features containing the different cues are Gaussian smoothed and interpolated onto the final simulation grid. The growth simulation uses continuous coordinates on a discrete grid.

The growing edges are implemented as instances of a growth cone class with the properties listed in Table 3. These include both internal parameters as well as functions that define how external cues are translated into growth directions. All object instances work on internal memories, are able to create new objects, and can inherit their parameters. The workflow and interaction of the different memories are summarized in Fig 4.

Table 3

Global parameters, attributes and actions of model agents.

Global parameter	Type	Description
Instances	i = 10	Number of global updates
Steps	s = 25	Number of growth steps between updates
MGD Switches	Boolean	Presence of different interactions
Agent property	Type	Description
Positions	[[x, y, z]i, ‥, [x, y, z]n]	History of all positions
Spherical angles	[[ϕ, θ]i, ‥, [ϕ, θ]n]	History of all directions
First and last node	list[str]	Initial node and node of last merge/branch/bifurcation
Current state	str	‘growing’, ‘merged’, ‘terminated’ or ‘out of bounds’
Action	Details
Growth	Construct MGD_ang(μ→,Σ), draw next direction
Branch	Create new daughter cone at an angle, inherit from mother cone
Bifurcate/Split	Create two new identical cones at an angle, delete old cone
Merge	Merge current with nearby edge, create new node
Terminate/Reactivate	Inactivate/reactivate cone, set status to ‘terminated’ or ‘growing’

Table 3: Internal attributes of nodes and edges, data types, and description.

Fig 4

Diagram of the simulation.

Starting conditions are provided as image volumes and as list of initial positions and directions. During each time step, the steps in the grey box are performed in parallel for each cone. The four elementary event types during network growth occur with their respective probabilities (A), depending on the local environment. Afterwards, the next position of each cone is determined from the spherical probability distribution computed from the local environment (B). Finally, the computation grid (left) and network dictionary (right), which are kept in shared memory, are updated by each parallel process, such that they can provide up-to-date conditions for all cones in the subsequent step. Finally, the network graph is constructed, its topological properties computed, and the network is visualized. These analysis steps can also be performed periodically during the simulation to monitor network development on-the-fly.

Prior to deciding on the growth direction according to Eq 8 angular, the cone checks if there is a merge, branch, bifurcation or termination event. For merging, the cone at position performs a collision detection with all edges E of an ensemble n ∈ Ω that are within its search field of radius R. Ω includes all edges of a cluster specified through its members that can interact during t growth steps. A new node is created at the ith positions of the jth edge j ∈ Ω that fulfils the condition to be in shortest distance to the cone position.

The cone terminates afterwards and the memory attributes of all participating objects are updated. Branching and bifurcation are executed via object functions. The events occur randomly with their respective probability. This probability and the parameters for branching angle and direction as well as the bifurcation plane are coupled to the external parameters of the spatial environment. The same type of coupling is implemented for termination and reactivation processes. If an agent branches or bifurcates, a new node is created at the splitting position. In the case of branching, a new cone object emerges from the new edge, and the mother branch updates its memory. For bifurcation, the old cone terminates, and two or more new cone objects are created and grow in new directions with an oblique angle to each other.

As the framework relies heavily on the communication of all objects, an efficient parallelization strategy is important. Already existing edges are not only stored in the graph dictionary, but also generate an imprint in an additional layer of the simulation grid. During growth, each process checks for the presence of existing edges by generating a view in this grid. If there is an edge in the vicinity, it can be looked up in the dictionary using its identifier from the simulation grid. To additionally reduce the number of collision calculations between actively moving cones, a Voronoi tessellation is performed for all growing edges to cluster neighbouring cones within a Delaunay distance l. Only those agents that are able to interact during the following l-steps have to communicate during these steps, while all clusters are independent from each other and can be updated in parallel. The clustering step is repeated after l steps. In Fig 5, a 2-dimensional clustering is shown as an example. All agents with overlapping blue circles can interact and are part of the same cluster.

Fig 5

Clustering and performance.

In (a), a 2-dimensional illustration of the clustering scheme to reduce the number of collision detections is shown for illustration. Blue dots correspond to active growth processes, and only processes with overlapping light-blue circles can interact. A Voronoi tessellation of the simulation volume is performed for all growing edges to cluster all neighbouring cones within a Delaunay distance l. Only agents that can interact during the following l-steps have to communicate during these steps, while all clusters are independent from each other and can be updated in parallel. In (b) and (c), the scaling of the simulation time with the number of growth cones and worker processes is shown.

Applications

We used our model to simulate the development of a multicellular biological network in a tissue with different degrees of anisotropy. This scenario is representative for a variety of biological systems where the final network architecture depends on the interplay of genetically encoded cell behaviour and the complex structured tissue environment. Examples include the osteocyte lacuno-canalicular system in bone [48], the development and regeneration of neuronal networks [2], or the reticular cell network in lymph nodes [15]. To show how our framework can be used to study the relationship between local growth rules, tissue structure, and global emergent network architecture, we performed growth simulations for different ranges of branching probability and angle, growth cone aperture, directional memory, as well as for different degrees of tissue anisotropy. To account for the stochastic nature of the growth process, we created a sample of ten random networks with identical starting parameters and calculated mean values. We then quantified the influence of the varying parameters on the resulting node number and average node degree, clustering coefficient, average shortest path length, edge density, as well as three types of node centralities. Closely spaced nodes are combined into a single node to enable comparison with experimentally obtained imaging data with finite spatial resolution [14]. Examples of the resulting networks and parameter relationships are shown in Fig 6, while the results of the sensitivity analysis of global network topology on the growth parameters are shown in Table 4. As can be seen in Fig 6A, the average node degree of the network gradually increases during the development of the network. A giant component emerges as all cells become connected to each other, indicating a critical point after around 30 steps. With increasing branching probability, the network becomes more dense, with higher mean degree and clustering coefficient and reduced average shortest path length (Fig 6B). Average node betweenness and information centrality decrease while harmonic centrality increases with higher branching rate. All three centralities approach the experimentally measured values in osteocyte networks [49]. We also find that stronger tissue anisotropy leads to a lower number of nodes but increasing clustering coefficient and average shortest path length, while the average node degree remains unchanged (Fig 6C). These observations suggest that the ability of growing osteocyte processes to branch at a sufficiently high rate together with a sufficient degree of tissue organization are important parameters for the formation of viable osteocyte networks in bone. This example outlines the possibility for more systematic studies, e.g. to identify the underlying molecular and cellular parameters for observed changes in network architecture in different biological systems. A similar comparison could be performed for any other 3D biological network for which such quantitative data are available. By modifying the provided code used to generate Fig 6 and Table 4, these examples can easily be adapted to other biological systems.

Fig 6

Example simulations of multicellular network growth.

a) Cells (blue) emanate growth processes that form a dense, interconnected network of active (green) and terminated (red) edges. Node degree goes up (top) and a giant component emerges that connects all cells (bottom), as evident in the connectivity graphs (red/black). b) Higher branching probability leads to increasing clustering coefficient and shorter average paths. Node betweenness and harmonic centrality approach experimental values of osteocyte networks in bone (blue). c) Stronger tissue anisotropy leads to smaller networks with increased clustering coefficient and shortest path length, and decreased information centrality, approaching experimental values (blue).

Table 4

Sensitivity analysis.

Parameter	N	〈k〉	〈cc〉	ASP	d	〈CB〉	〈CI〉	〈CH〉
Directional Memory	↗	-	↘	↘	↘	↘	↗	↗
Growth Cone Aperture	-	-	↗	-	↗	↗	↘	↘
Branching Probability	↗	↗	↗	↘	↘	↘	↘	↗
Branching Angle	-	↗	↗	-	-	↗	-	↘
Tissue Anisotropy	↘	-	↗	↗	↗	↗	↘	↘

Table 4: Sensitivity analysis of global network topology on local growth parameters and the anisotropy of the tissue environment. Columns from left to right: Number of nodes, average node degree, average clustering coefficient, average shortest path length, edge density, node betweenness, information and harmonic centrality.

Discussion

We developed a novel agent-based framework to simulate the growth of networks in three-dimensional space. The outgrowth of individual processes follows a 3D random walk with memory (correlation) and directional bias due to interactions with the environment. This framework can be used to simulate the growth of biological networks on all scales from single cells to ecosystems, or of generic network development through interaction of autonomous agents with each other and with external physical and chemical cues. During simulation, the network is directly converted into a graph using the networkx library [50], which provides access to many quantitative and comparative measures of network topology.

The first agent-based computer simulation of 2D spatial network growth, published in 1967 [20], already used a similar approach as our model: edges grow and branch according to rules, respond to attractive and repulsive cues, and follow density gradients. The next growth direction is randomly assigned based on sampling the space around the growth tip. Inspired by plant growth, this model can generate a variety of tree-like networks—a remarkable achievement given the limitations of computer simulations more than 50 years ago. While the field of complex networks flourished decades later, relatively little work was done on spatial network growth [1, 19]. In 2009, two new computational frameworks for neuronal network growth were published, CX3D [25] and NETMORPH [26]. Both use biologically realistic growing and branching rules and are openly available. In both models, synapse formation only takes place after the growth phase. CX3D offers a more physically realistic treatment of mechanical forces compared to our model, but lacks the flexibility of grid-based external structural cues. Two other neuron growth models [23, 51] efficiently generate neuron tree morphologies and space-filling networks, but without biologically realistic growth mechanisms and synapse formation. Taylor-King et al. [52] use an elegant mean field approach for network growth based on local state degree distributions. Whenever it is biologically plausible to analytically formulate local rules, such an approach can accurately reproduce global properties. In their model for vessel tree development, Perfahl et al. [53] include physical interactions with the environment and between edges to study the role of mechanical forces for vessel development. In their “Unifying theory of branching morphogenesis”, Hannezo et al. [28] present a stochastic model for epithelial duct development based on branching and annihilating random walks. Although interactions are limited to explicitly defined chemical gradients and anisotropy, this model can reproduce experimental tree topologies for a variety of organs. We summarized this comparison of our model with prior work in S1 Table. To our knowledge, our work presents the first framework for network growth that includes a biologically motivated local growth process as well as arbitrary interactions with an external environment, and is fully available as computational implementation.

We did not include energy consumption in our model in order to keep it as generic as possible. Depending on the type of network, different mechanisms of energy sources, energy distribution and consumption are possible. As an example, energy could either be “harvested” from the environment by the growing processes, or distributed throughout the network from cell bodies. Other existing models either omit energy consumption as well, or contain only specific energy-related mechanisms such as competition between growth cones [26] or diffusion of proteins from the soma [25].

One class of biological networks where the interaction and feedback between invidiual growth cones and the microenvironment during growth determines the resulting network architecture are multicellular networks, such as neuronal networks in the brain or the osteocyte network in bone. In bone, different degrees of tissue organization and the presence of soluble cues determines the connectivity and arrangement of the resulting network. With the framework presented here, the question can be investigated how local cellular response and resulting network topology depend on each other by conducting virtual experiments with varying parameters, which would not be possible in real experiments. The resulting theoretical understanding can then inform the design of guiding structures (scaffolds) to facilitate optimal network development in regenerating bone by providing the relevant physical and chemical signals [54].

The growth of vascular and neuronal cell networks during embryonal development is another example where soluble cues, the geometry of the environment, and the interaction between nearby cells all play an important role [55]. Recent improvements of lightsheet imaging technology now allow to visualize the entire growth process with high spatial and temporal resolution [56]. Our framework predicts not only the outcome, but also the dynamics of this process, which makes it possible to test specific hypotheses about the biological mechanisms that locally control this growth process. In the future, this might contribute to solving the question to what degree functional network architecture is encoded in the genome, or how the interactions and mechanochemical feedback loops between cells and the environment during growth determine the resulting tissue organization. Another interesting example is the role of curvature for tissue growth and organization [57]. With our framework, it will be possible to systematically vary the curvature of growth surfaces for networks while including other relevant biological and physical signals in the same simulation framework.

A very interesting related field is the exploratory growth of plants guided by pairwise interactions, tropic reponses to signals (e.g. light), and nastic (e.g. helical) movements. Recently, a computational model was introduced to model the dynamics of such sensory-growth systems [58, 59], taking into account the mechanical properties of the growing system. Such a framework offers interesting opportunities not only to understand biological control principles, but also for designing self-growing artificial systems. It would be interesting to see how this approach can be extended towards the development of functional network architectures.

While we developed our framework with network growth in mind, it can also be turned around and used to find and trace network structures in noisy image volumes by treating them as guiding signals for growth. Many approaches exist to trace filaments [60, 61], but they often have difficulties with branching structures. By simulating multiple growth processes either directly on image data or e.g. on probability maps predicted by convolutional neural networks, our framework might be able to find also highly branching and irregular structures.

Network growth as defined here is a parallel process by definition and thus can be computed in a distributed manner. We describe an efficient strategy for parallelization while taking into account the interaction between neighbouring growth cones. Future improvements could make use of the massively parallel processing power of graphical processing units.

Conclusion

We developed a probabilistic agent-based model to describe individual growth processes as biased, correlated random motion with rules for branching, bifurcation, merging and termination. Using directional statistics, the influence of external cues on the growth process could be obtained directly from real image data. Our tool allows us to study the relationship between biological growth parameters and macroscopic function, and can generate plausible network graphs that take local feedback into account as basis for comparative studies using graph-based methods. Our approach is not limited to a specific type of network, includes a thorough treatment of probabilities, and can easily include arbitrary external constraints. We implemented our model in python and make it freely available as a computational framework, that allows other researchers to perform simulations for diverse types of biological networks on all scales.

Methods

All edges and nodes are instances of their respective classes, with the attributes shown in Table 3. The framework was implemented in Python 3 using numpy [62, 63]. The network is directly converted into a networkx graph for further processing and analysis. On a 8-core Intel Xeon 6134 3.2 GHz machine with 64 GB RAM, running a network simulation with 100 growth cones in a 256x256x256 tissue volume for 50 steps takes 60 seconds. The scaling of the simulation time with the number of growth cones and worker processes is shown in Fig 5. Source code, documentation and example configurations of our computational model are available at https://github.com/CIA-CCTB/pythrahyper_net.

Detailed description of 3D random walks in complex environments.

(PDF)

Click here for additional data file.

Extended version of Table 1 with different external signals and cues.

(PDF)

Click here for additional data file.

Comparison of our framework with prior work.

(PDF)

Click here for additional data file.

27 Jul 2020

Dear Prof. Kollmannsberger,

Thank you very much for submitting your manuscript "Biological network growth in complex environments - a computational framework" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Please note the attached reproducibility report, which should help you improve the reproducibility of your manuscript.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Pedro Mendes, PhD

Associate Editor

PLOS Computational Biology

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Reproducibility Report has been uploaded as an attachment.

Reviewer #2: This manuscript addresses a major gap in our understanding of complex life: exploring generative processes which give rise to multicellular assemblies. It is very relevant work to an under-researched area.

An agent-based model is presented which gives rise to complex branching structures. The model is well described. The outputs of this generative process appear to be interesting and potentially relevant.

There are several areas which can be expanded upon to strengthen this manuscript:

- The biological relevance of the generative model could be made more explicit. In particular, how the various parameters in the model might be implemented or represented in a biological system. This would help support the inclusion of these parameters and being tangible modulators, rather than a series of variables capable of reproducing biological observations, but detached from how these processes actually occur.

- The networks generated by the model visually appear similar to the 2 biological networks they are compared to in the study. However, no quantitative comparison is provided. Using a variety of centralities to quantitatively compare the generative models to biological data would strengthen the conclusion that these models are capable of creating biologically meaningful outputs

- Following on from the previous point, a sensitivity analysis as to how the various parameters in the model impact the topological properties of the outputs from the model would be very useful. Again, using quantitative approaches via centrality measures, the range of topologies generated by the model, and the presence of critical points in parameters can be identified.

Reviewer #3: Major Comments

I enjoyed reading this paper, and I believe the premise of the paper - a generalized model for biological network growth under a variety of conditions - is extremely promising. The authors put great thought and effort into deriving the mathematical principles of network growth (in particular, I thought it was very clever to model individual environmental cues as multivariate Gaussians that can be combined through convolution into a single aggregate "environment" constraint). I am by no means an expert in directional statistics or geometry (and I have some points of confusion described below) but from what I understood, the mathematical foundations of their network growth model appear sound.

While the methods presented in this paper are exciting, I have concerns about publishing this without more empirical results. The authors present two anecdotal examples of using their network growth model to generate networks similar to those observed in biology. These two examples are not very convincing; in Figure 6, it is not obvious to me that the simulated network actually replicated the original network. And even if it did, these are just two examples. The authors are advertising a mathematical framework with a robust ability to generate a wide variety of biological networks by simply tweaking the input parameters, but they don't provide very much evidence of this. I realize this is easier said than done, but I think the paper would be strengthened with quantitative evidence that their network growth model can reproduce the growth seen in a variety of systems. This would likely involve finding a particular set of parameters that, when used in simulations, generates a sample of random networks that are then shown to be quantitatively similar to observed networks from the same system (and then repeating this process for other systems).

Minor Comments

- typo: "our tool allows to study..."

should this be "allows us to study..."? Note that this phrasing occurs throughout the paper, and if it is changed here it should be changed throughout.

- line 51-52: "making too many assumptions about these cues"

- line 270-271: "not modelled explicitly, e.g. in terms of hard-coded conditions, but through a generic mechanism based on angular probability density functions."

I think the authors would do well to more explicitly explain how their model for directional

cues differs from prior work.

- are equations (1) and (2) necessary? As the authors state, an analytic solution is elusive, and the authors end up resorting to a numerical solution based on random walks.

- I'm having some trouble understanding lines 82-90. In particular, how should I interpret the "mean vector" \\mu? Does this encode the location of some signal? Or does it encode the direction in which the current edge is being "pulled"? If it's the former, I think this needs to be better explained. If it's the latter, I believe it would be helpful to include a subscript, i.e. \\mu_t, to indicate that this is a time-dependent value.

This confusion carries on throughout the paper, and I think it could use a more intuitive explanation in lines 82-90.

- line 96 - should "determinate" be "determinant"?

- why does the document intermittently stop listing line numbers?

- I am slightly confused about the section on projecting onto the unit sphere. My intuition is that we essentially generate a random direction vector, and then project it onto a sphere of a randomly chosen radius? Is this correct? If so, I think the paper would benefit from explaining this explicitly.

- the MGD in equation (3) takes three inputs; but the MGD inside the integral in equation (7) takes only two inputs; why the disparity?

- line 162: should there be seed directions, in addition to seed locations?

- I could not understand the purpose of representing the environment as a grid. Is this simply for clustering to efficiently compute collision and merge events? Otherwise, I am not sure what difference this makes since the coordinates of edges are still continuous.

- Figure 6: are the blue/green square drawing, and the blue circular drawing, both graphical representations of the same simulation? This was confusing to me

Reviewer #4: The author developed a computational model of biological network growth in complex environment. A computational framework is based on directional statistics to model network formation in space and time under arbitrary spatial constraints. Growth is described as a biased correlated random walk where direction and branching depend on the local environmental conditions and constraints. To demonstrate the application of the computational model, authors performed growth simulations of the osteocyte lacuno-canalicular system in bone and of the zebrafish sensory nervous system.

It does a nice job of introducing relevance and background for such a question and overall I think it is interesting computational material. Unfortunately, a couple of things have put me a bit off enjoying the quality and impact of the paper's results:

1. The structure of extracellular space is named “images”. What is description of “images”? A more detailed description of the «image» concept must be given.

2. The variable p is introduced in the equation 1, but it is not written what this variable means in the framework of the model.

3. On page 5 (line 78) the vector u is indicated in bold, although throughout the entire article vectors are indicated by an arrow at the top. It is necessary to use the same designations throughout the article.

4. It is necessary to describe a possible biophysical mechanisms of “self-avoidance”, and references to experimental works must be provided.

5. What does means parameter “memory” in Table 2? There is no specific description or explanation in the text. Does this mean computational memory?

6. The model does not take into account energy consumption of growing branch at all. The authors should to explain in the text of the article, why this important process (which can create restrictions) does not included into the model.

7. For what purpose is table 4 in the article? This is purely a description of Python data types. I think this table needs to be dropped.

8. Authors used the model to simulate the development of two biological network systems: the osteocyte lacuno-canalicular system in bone and the sensory nervous system in Zebrafish. But they obtained only qualitative conclusions that states, that there is growth, and the images obtained are similar to experimental pictures. I believe that it is necessary to provide some quantitative characteristics that would allow readers to estimate how close the simulation results are to experimental results (for example, maybe by using the metrics of random graphs theory).

9. There are various models that describe growth in biological systems and networks. Authors should provide a brief comparison of their model with existing models, and should specifically indicate the advantages of their model over others.

10. In «Conclusions» (line 327) a reference to a specific paper(s) must be given.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: None

Reviewer #2: Yes

Reviewer #3: Yes

Reviewer #4: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Anand K. Rampadarath

Reviewer #2: Yes: George Bassel

Reviewer #3: Yes: Arjun Chandrasekhar

Reviewer #4: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

Submitted filename: Reproducible_report_PCOMPBIOL_D_20_00894.pdf

Click here for additional data file.

2 Oct 2020

Submitted filename: Response_to_Reviewers.pdf

Click here for additional data file.

29 Oct 2020

Dear Prof. Kollmannsberger,

We are pleased to inform you that your manuscript 'Biological network growth in complex environments - a computational framework' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Pedro Mendes, PhD

Associate Editor

PLOS Computational Biology

Jason Papin

Editor-in-Chief

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #2: The authors have addressed all concerns raised

Reviewer #3: No further comments

Reviewer #4: The authors have done a serious job. The article is now much better than the first version. I believe that the article can be accepted for publication.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #2: Yes

Reviewer #3: Yes

Reviewer #4: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: No

Reviewer #3: Yes: Arjun Chandrasekhar

Reviewer #4: No

24 Nov 2020

PCOMPBIOL-D-20-00894R1

Biological network growth in complex environments - a computational framework

Dear Dr Kollmannsberger,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Nicola Davies

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol