Critical Issues in Modelling Lymph Node Physiology

In this study, we discuss critical issues in modelling the structure and function of lymph nodes (LNs), with emphasis on how LN physiology is related to its multi-scale structural organization. In addition to macroscopic domains such as B-cell follicles and the T cell zone, there are vascular networks which play a key role in the delivery of information to the inner parts of the LN, i.e., the conduit and blood microvascular networks. We propose object-oriented computational algorithms to model the 3D geometry of the fibroblastic reticular cell (FRC) network and the microvasculature. Assuming that a conduit cylinder is densely packed with collagen fibers, the computational flow study predicted that the diffusion should be a dominating process in mass transport than convective flow. The geometry models are used to analyze the lymph flow properties through the conduit network in unperturbedand damaged states of the LN. The analysis predicts that elimination of up to 60%–90% of edges is required to stop the lymph flux. This result suggests a high degree of functional robustness of the network.


Introduction
The immune system functions to protect the host against various pathogens and tumours. The mammalian immune system is a complex compartmentalized population of different types of humoral factors and cells organized as a body-wide network of interacting lymphoid organs, e.g., bone marrow, thymus, spleen, lymph nodes, gut. Its elements are connected by two vascular systems, i.e., the blood and the lymphatic systems, that enable the immune elements to migrate continuously to or from the compartments during their lifetime. Lymph nodes (LNs) represent the most numerous compartments of the immune system, in which soluble signals and cells are concentrated to enable efficient interaction between them. The structure of a LN is crucial to its function [1]. The development of anatomically correct mathematical models that accurately reproduce the physiological, cellular and molecular processes in LNs remain one of the major challenges for mathematical immunology.
Modelling the lymphatic system and in particular, the lymph flow through LNs is in its infancy [2]. There exists very few computational studies in which the fluid flow through the internal domains of LN was investigated [3,6]. The model is based on an idealized 3D geometry of the LN macro-structures and an approximation of the B-cell follicles, T cell paracortex and medulla as porous 2 of 22 regions. The multi-physics description of the lymph flow processes combines the conservation of mass and momentum equations in subcapsular and medullary sinuses (fluidic region) with the conservation of mass and Darcy's law with Brinkman's term for the flow in porous regions. This study can be considered as a milestone for the application of computational fluid dynamics and the whole-organ fluorescent imaging technologies to model the LN physiology.
Lymph nodes are characterized by a highly variable histological appearance. This implies that the development of computational models of the LN structure and function requires (i) an idealization of the structure, (ii) the implementation of a modular approach, and (iii) the analysis of the performance of specific functional LN units in health and disease. The elements of the above framework can be found in [3][4][5][6][7].
LNs are known to present high resistance to lymph flow [3], which results from the internal structure of the LN body comprised of blood microvascular networks, the conduit networks and of domains densely packed with immune cells, i.e., the B cell follicles and T cell zone. In the existing 3D geometry LN models developed for studying the fluid flow, the conduit and blood vascular systems are not structurally described. Meantime, these elements, in particular the conduit system network ensheathed by fibroblastic reticular cell (reticular) network, represent a key system for the delivery of soluble molecules in the paracortex area of the LN. It was estimated in [3] that under normal conditions about 10 % of lymph entering the LN through the afferent vessels takes a central path to reach the efferent vessel and high endothelial venules (HEV) via the conduit system and the inner domain of the LN. Under inflammatory conditions and during infections, the conductivity of the conduit system is likely to change dramatically [8], which affects the information delivery (cytokines, hormones, etc.) to the reactive parts of the LN. A fine and elaborate structure of the LN vascular systems requires the development of multi-scale computational models to comprehensively examine the fluid transport in LNs.
In this study we present an overview of the critical issues in modelling the structure and function of LNs and propose some approaches to model the networks of fibroblastic reticular cells, the blood microvasculature and analyze the flow properties of lymph through the conduit network in an unperturbed state and a damaged state with the varying degrees of destruction. We review the major structural elements of LNs in Section 2. The approaches to model the geometry of an a FRC network and the blood microvascular network are presented in Section 3. A fine resolution analysis of the lymph transport through a single conduit are presented in Section 4. The hydraulic conductivity properties of normal and disrupted conduit systems, modelled as 1D topologically equivalent networks are studied in Section 5. Numerical assessment of the percolation robustness of the fibroblastic reticular cell (FRC) network is predicted in Section 6. Further perspectives for application of computational fluid dynamics techniques are discussed in Section 7.

Major Structural Elements of a Paradigmatic Lymph Node
Lymph nodes (LNs) are specialized secondary lymphoid organs which are strategically situated at sites where pathogen invasion is likely to occur and where immune responses are induced and maintained [1]. The idealized LN structure can be divided in three major regions: the outer cortex containing the lymph-draining subcapsular sinus (SCS) and B cell follicles, the inner paracortex which constitutes the T cell zone and the medulla containing the efferent lymphatics. Each region contains specialized cellular niches which are critical for the development of adaptive immune responses [9]. B cell follicles are located underneath the SCS where B-cell specific responses take place [10], whilst the inner LN area constitutes the T cell zone where T cells are activated through interaction with antigen-presenting cells, such as dendritic cells (DCs) [11]. Fibroblastic reticular cells (FRCs) form a densely intertwined network that supports T-DC interactions, subsequently leading to generation of immune responses [12][13][14]. Furthermore, FRCs support lymphocyte migration, retention and survival through expression of chemokines CCL19 and CCL21 [15] and survival factor IL-7 [16,17]. FRCs secrete a basement membrane enwrapping the inner core of LN conduits which consist of collagen fibers and extracellular matrix (ECM) components. The conduit system drains lymph from the SCS through the entire LN T cell zone and the interstitial lymph flow accommodates small molecules and antigen that are < 80 kDa in size. Some DCs are capable of extending their cell protrusions into the conduits, enabling them to sample the inner core for foreign antigens [18]. Apart from the lymphatic system in LNs, a specialized blood vascular system permeates through the LN T cell zone, forming high endothelial venules (HEVs) [8,19] which serve as extravasation points for migrating T, B cells and DCs into the LN [20]. Despite technical limitations in terms of rendering difficulties in the assessment of structural information from whole LNs, recent advances in imaging and computational techniques have enabled a comprehensive quantification of e.g. the LN blood vascular network [8,21].
The complexity of the immune system has been recognized in the recent decade as one of its key hallmarks, but also a major challenge for integrative immunological research [22]. Nevertheless, the development of high throughput "-omics" methods [23] and systems biology approaches [24] has made it possible to dissect the spatiotemporal interactions between immune system constituents and reveal functional implications in the context of various diseases. In order to fully utilize systems approaches, it is necessary to integrate available data across multiple biological scales, from lower level processes such as gene expression in different cell populations to higher level processes such as maintenance and regulation of tissue and organ functionality [25]. Finally, the complex dynamics of the immune system can be elaborated from the content-rich observations of these processes.
Increasing computer processing capabilities has made it feasible to design and generate a paradigmatic LN model based on real organ geometry [4]. Such a systems biology approach has allowed us to develop a generalized 3D solid model of the LN, which consists of the following idealized structures: SCS, B cell follicles, trabecular and medullar lymphatic sinuses, blood vasculature and HEVs, the T cell zone and FRC network. An integrated LN model across multiple scales will allow detailed description of global spatiotemporal dynamics of the different cell type interactions and enable us to identify biological parameters that are critical for generation and maintenance of immune responses.

Cellular Potts Modelling of the FRC network
In this section we construct a geometrical model of reticular network using the Cellular Potts Model (CPM). This approach lets us explicitly control the target volume of reticular network, vary properties of each FRC, and directly incorporate the network in further models of cellular immune response. In addition, it is a starting point of modelling reticular network as soft tissue, which would let us explore mechanical properties and behavior of the network in cases of inflammation, T cells depletion or LN fibrosis.
In CPM (reviewed in [26][27][28][29]), the cells are defined on a cubic lattice as connected collections of lattice cites (voxels) v = {i, j, k} with the same cell id σ(v). Each cell can belong to a certain type τ(σ(v)). The extracellular medium is defined as a generalized cell with σ = τ = 0.
The dynamics of the system are defined with a modified Metropolis algorithm, which models the motion of cells by rearranging the states of the voxels in a stochastic energy minimization manner. It samples random local voxel-copy attempts and accepts them in accordance with provided Boltzman probability acception function p(σ(v source ) → σ(v target )).

of 22
The energy of the system is usually defined as the sum of two terms: Adhesion term represents the sum of binding energies J(τ(σ(v 1 )), τ(σ(v 2 ))) of molecules located at adjacent pairs of voxels v 1 , v 2 , which correspond to membranes of different cells σ(v 1 ), σ(v 2 ). The second term describes the ability of the cell to constrain its volume V(σ) near the resting volume V target (σ) due to the internal pressure. The parameter λ(σ) is a spring modulus, i.e. the rigidity of the cell to changes of volume. It determines the relative weight of the second energetic term. There are many other possible terms and modifications of the form of energy, E, which allow the modelling of compressibility of cells membranes, chemotaxis, haptotaxis, cell elongation and other processes [26,27,29,31]. The basic form of the probability acceptance function in CPM is given as, where T is a global parameter representing the amplitude of cells fluctuations, proportional to the probability to accept energetically unfavourable voxel-copy attempts, ∆E = E a f ter − E be f ore is the change of energy of the system (1) as a result of sampled voxel-copy attempt. If T = 0, the simulation is fully determinate and can freeze at local energy minima. As a possible extension of (2), each cell type or each cell can have its own T(σ).
The other form of acceptance function, which provides a biologically relevant meaning of parameter T(σ) -the intrinsic motility of the cell (discussed in detail in [29]), was specified as follows, where is the resultant motility of membrane between σ(v source ) = σ(v target ). Given the tanh-mapping of membrane motility, function (3) restricts the voxel-copy attempts involving so-called frozen cells (with T(σ) → 0), even if calculated ∆E is significantly negative. This allows us to specify the motility of cells, i.e. their ability to react on external stimuli. We extend this approach by setting intrinsic motility for each voxel of the cell. By doing this, we can easily model the reticular network, making the ends of protrusions frozen (forming FRC junctional complexes [30]). As a result, there is no need for extra assumptions and extra CPM plugins such as cell elongation or connectivity constraints (which are computationally expensive in 3D [31]), utilized in other models of cellular networks (such as in the study of vasculogenesis [31]).
We initialize cells as thin (d ∼ 1.0 µm) cylinders along the graph of reticular network topology (constructed in [4]), forming the conduit system ( Figure 1a). Protrusions of different cells are connected at the middles of the corresponding edges. For each FRC we define the center of it's body v c (σ) at the coordinates of corresponding graph nodes. We define the motility of voxels occupied with cells (v(σ)) according to the distribution localized around v c (σ): By doing this, the voxels, which are far from the center of FRC body v c (σ), have close to zero motility T(v(σ)) → 0. As a result, the ends of FRCs protrusions are linked with each other in accordance with given topology and are not altered during the simulation, resulting in no connectivity issues. The voxels which are closer to centers v c (σ) have more intrinsic motility, thus it is more probable for them to extend to the medium due to internal pressure until the FRCs would reach their target volume V target (σ). Figure 1b illustrates the reticular network, reached it's target volume (4% of the whole volume of lattice [4]). It contains 3373 FRCs in a 177.0 µm × 197.4 µm × 201.0 µm computational domain with l px = 0.3 µm length of the voxel resolution. Parameters used in the simulation are listed in Table 1. We set adhesion energies to zero because they are not relevant for the result as there are no other cells in our model. The choice of FRC-FRC adhesion energy doesn't affect the simulation because FRC junctional complexes are modeled by freezing the ends of protrusions. A high value for lamda is chosen (λ(σ) 1) because of rigidity of the reticular networks. The deviation of volume from the target volume for the whole network is less than 1% after 20 Monte-Carlo steps (MCSs) burn-in period.  The implementation of the model is based on the open-source CompuCell3D application (available at http://compucell3d.org/). We modified the source code of its core C++ library to accomplish voxel-based motility and developed Python scripts to configure the simulation. The computation of the FRC network shown in Figure 1 requires about 10 min real time on Intel Xeon 4-core 3 GHz CPU with 4 Threads based parallelization.

Modelling Blood Vascular Network
The vascular network of a lymph node is a rather complex construct and is placed inside the internal space of the LN with various items, such as B-cell follicles, FRC network, medullar and trabecular zones. The vascular network provides the delivery of oxygen, nutrients and signal molecules. Lymphatic fluid can leave the LN through the blood vessels due to a pressure difference, hydrostatic and osmotic, between the interstitial pressure in the node and the pressure in the blood vessels [2,3]. The vessels of the blood microvascular network have lots of close intersections with the FRC network, pervading the whole space inside LN, so the vascular network influencesshma many processes in the lymph node. In this way, construction of a 3D model of vascular network provides an insight into the role of the vascular network in processes of tissue homeostasis and pathology.
While constructing the vascular network, we had to solve some problems with intersections of the FRC network, vascular network and internal items of the an LN. The main challenge there was to overcome the differences between the typical sizes of the networks. An FRC network consists of conduits with a diameter range from 0.2 to 2 µm, while the vascular network contains vessels with a diameter range from 5 to 30 µm [21]. We modified the algorithm from [7] to meet these conditions.

Initial Data
We used data from [21] to set the length of the network vessels. The blood microvascular system  Blood vessels length distribution summarized from [21] is characterized by the length of vessels and branch separation from the feeding vessel. The number of capillary vessels was taken from [21] as 40 vessels. It could be approximated as 3 vessels with 5 bifurcations (32 vessels) and 1 vessel with 6 bifurcations (64 vessels), the diameter of input and output vessels was taken as 8 µm.

Algorithm of Network Graph Generation
The algorithm of vessels network graph construction consists of three main steps: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 9 December 2016 doi:10.20944/preprints201612.0054.v1 Step 1. Graph topology organisation. In this step we generate the basis points and edges of connections.
Step 2. Local edge length optimization. In this step we use the algorithm from [7] just for a local (i.e., for neighbouring nodes) adjustment of the mismatch of the model and target graph edges lengths . In this and the next steps, the following parameter from Step 1 is used: blos (defined on line 4 in the pseudocode). It's the canonical length of segments of the vessels graph.
Step 3. Global network structure optimization. In this step we use a modified algorithm from [7] for (i) minimization of the edge length deviation from the real data for all neighbouring nodes, (ii) pushing apart disconnected nodes from each other to prevent merger of the vessels, and (iii) shifting the nodes away from the prohibited domains associated with other LN structures.
Pseudocode of Step 1 (graph topology organisation): // Initialise the data arrays 1 set "length distribution" as real array ld; 2 set {5 5 5 6} as integer array bf; // Specify the segmentation accuracy, vessels radius and decreasing, processing zone size 3 set "initial input/output vessels diameter" as real vd; 4 set 4.0 as real blos; // length of segmentation of the vessels, µm 5 set 2 1 2 as real rf; // coefficient of vessels radius decreasing 6 set "work sphere radius" as R; // simplify the constructing 7 set graph structure t; // Attach the input and output vessels 8 insert line "input vessel" and line "output vessel" to t; // simple lines, splitted into 100 segments 9 f or (integer j = 1 to NC) begin // In this loop, we create new vessels, growing from input and output vessels 10 set real vdl = vd; 11 set real sl = random f rom ld;    The blood microvascular network generated using the developed algorithm is shown in Figure 3.

Integrative Geometric Model of Vascular Networks
To assemble the two networks [FRC + vascular networks] into one 3D geometric model, we used the FRC network constructed with the algorithm from [7] for a sphere with a diameter of about 200 µm as shown in Figure 4. The graph of the network topology and its node-level characteristics are shown in Figure 5.  The method developed in [7] was further utilized to place two or more network graphs in same domain. To this end, we placed both graph structures in one structure, unified the length of all edges (we converted them into sets of edges with length about 4 µm) and executed the minimisation of length inconsistencies using a modified version of the algorithm from [7]. Some remaining intersections were removed at the stage of voxel-based approximation of the graphs. The key characteristics of the constructed 3D blood microvascular model are listed in Table 2:

Lymph Dynamics in Conduit Elements of FRC Network
In this section we examine generic conduit and the FRC network properties.

Transport Through a Single Conduit
FRC network conduits consist of reticular fibres arranged in bundles encased by FRCs [30]. These conduits have been observed with radii between 200nm and 3µm [35], consisting of 10-100 reticular fibres [36] with radii of 20-40nm [37]. Microscopy of FRC conduits appear to show that the fibres are densely packed [36]. Lagrange found the optimal packing of cylinders is hexagonal with a maximum fibre density of φ = π 2 √ 3 . There is some justification for the assumption that FRC conduits consist of reticular fibres arranged in such a manner. As the FRCs hold the fibres together it is likely that they exert a force on the fibres pushing them close to one another. The above numbers would suggest that a 200nm radius conduit would have between 7 and 29 fibres which would appear to approximately agree with values found in literature [36]. A diagram showing the paths between fibres available for flow and the arrangement of the fibres is shown below in Fig. 7.
The Navier-Stokes equation for incompressible, steady-state, unidirectional flow of a Newtonian fluid reduces to a Poisson equation and can be written as below.
Where u is the velocity, µ the viscosity and p z is the z axis component of the pressure gradient. Using a finite element solver with a no-slip boundary condition u(Γ) = 0. An area independent flow parameter, χ, is defined as follows, A finite element solver was used the find χ for curvilinear triangular domains such as Ω. The area independence of χ was verified and a convergence study performed, the results are tabulated below. A rearrangement of eq. 7 gives the flow in one of the Ω domains for a given pressure gradient in the form QR = p z with R defined as , Based on the values found for χ it is observed that the resistance value found for flow through a curvilinear triangular cross-section is twice that through a circular cross-section of the same area. The total resistance, R tot , for a conduit containing multiple fibrils is the reciprocal of the sum of the reciprocal individual resistances provided that r c >> r f . This leads to R tot = R/n, where n is the number of fibres.

Diffusive Transport in an FRC conduit
At scales such as these diffusive mass transport can be considerably more significant than fluid flow. The geometry of the channels seems that it will produce large viscous losses whilst still retaining a large cross-sectional area available for transport. Assuming ideal mixtures a relationship (Fick's law) can be constructed for mass diffusion within the conduits in a similar manner, where ∆C is the difference in molar concentration, N is the molar flux and D is the diffusivity. Now if we assume a 300µm long pipe spanning a lymph node with a pressure difference of 6 cmH 2 O and a viscosity of 1.5cp [38] then the Péclet number in a conduit is 0.0388 -for diphtheria toxin, see below -suggesting its motion is diffusion dominated. Thus it is reasonable to assume the flow will have an insignificant effect on the concentration of diphtheria toxin allowing these two transport phenomena to be decoupled. Instead of a 3D fluid solver operating on a void space of the network, a 0D solver can operate on the connectivity and associated properties of the network with a considerable reduction in computational expense. Such a solver was implemented in MATLAB R and a cuboid block of a generated FRC network consisting of 6927 edges and 3374 nodes was subjected to analysis. Contiguous nodes were made co-planar to give a domain of known dimensions, 151 × 193 × 187µm. Pressure or concentration gradients were placed across the block, a linear system formed and solved for the fluid permeability or diffusive flux using a direct solver. One species often considered in lymphatics is diphtheria toxin with a diffusivity D = 6.2 * 10 −7 cm/s [39] and a molecular weight of 58kDa [40]. For the network under consideration a 2 ng/g imposed gradient across the first axis of the network gave a diffusive flux of 2.68 µg/s. Whilst for an imposed pressure gradient of 6 cmH 2 O gives a fluid flux of 0.0214 µm 3 /s. The fluid flux would have to be ∼ 10 orders higher to have a similar contribution to mass transport as diffusion.

Normal FRC network
In order to simulate a steady-state lymph flow in a conduit system of the FRC network model in an idealized LN designed in [7], we applied Poiseuille law as well as mass conservation law. Let's consider the graph of conduit network shown in Figure 4 and denote set of vertices and edges G = (V, E). For each edge e ij ∈ E we apply the Poiseuille equation and for each vertex v i ∈ V we apply the mass conservation law (i, j = 1, 3694).
The variables are the lymph flow Q ij [(µm) 3 /s], the hydraulic pressure P i [Pa] and hydraulic resistance, R ij .
According the Poiseuille equation, hydraulic resistance for non-elastic tubes is the following: R = 8µl πr 4 The parameter µ = 0.0015 Pa·s is the lymph dynamic viscosity [32]. The radii r ij of all channels are assumed to have a constant value of 1µm. Channel lengths l ij can be calculated from graph data shown in Figure  4. If we take into consideration reticular fibres densely packed in conduits and use hydraulic resistance defined in equation 8, the flow through the system will be about 7 orders lower.
The system 10 needs to be closed by setting boundary conditions. The conduit network is connected

Gradient Constant
Pmax Pmin x Af Ef . Figure 9. Schematic view of the SCS pressure distribution for two types of pressure boundary conditions ((1) gradient and (2) constant. Branches emanating from the inner SCS circle represent the input and output edges of the network to the floor of the subcapsular sinus. In our model implementation, the nodes located at the border of the FRC network graph are connected to the SCS by additional edges. The flow through the system depends on pressures in these vessels. There is no experimental data on the pressure distribution in the SCS. In this study we consider two variants of pressure distribution (see Figure 9). If we set the x-axis in the direction of flow through the LN then boundary vertices with smaller x-coordinate will have higher pressure values. The first variant of pressure distribution is called «Gradient» with a linear decline from P max = 500 Pa to P min = 0 Pa depending on x-coordinate. The pressure estimates are based on values reported in [3,6,33]. The second variant is referred by term «Constant». If we set a plane normal to the x-axis dividing the LN into two halves of equal length, points in the upper half will have a value of P max and points in the lower half will have P min .
Having specified the boundary conditions, we build a linear system, excluding the flow variables from equation 10, and solve it for pressure. The system matrix is sparse but nonsingular so we use UMFPACK solver[34] which provides a direct method for linear sparse systems. The following results have been calculated using the UMFPACK direct method. The inflow and outflow were tested to be equal, the conservation of mass condition is satisfied. A network consisted of 3694 vertices and 7253 edges with the number of input and output vertices 169 to 151 in the «gradient» case, and 164 to 156 in the «constant» case.  Figure 10 demonstrates the change of pressure with growth of distance from the afferent lymphatic vessel for all vertices with different boundary conditions. In case of the gradient boundary condition (10a) there is an evident linear gradient with a line of boundary points from which others deviate slightly. In case of the constant boundary condition the pressure gradient is obviously higher, but still the pattern is close to linear.

Disrupted FRC network
Viral infections suffered by the immune system can destroy parts of the FRC network and lead to its disruption [8,14]. We applied the model of lymph flow in the conduit system to parameterize the destruction process and study its effect on lymphatic system function. To simulate a disrupted conduit system we delete edges from E randomly. If all edges connected to a node are deleted it becomes isolated, with zero pressure. Deletion of edges also leads to broken chains and pendant (or leaf) vertices. Mass conservation law implies that there is no flow into pendant vertex and no flow through broken chains. Some edges which are not removed may have no flow and so disfunction, we call these edges non-functional («n-f edges»). For example, in Figure 11 illustrating the disruption elements vertices 2 and 3 are pendant, edges (1-2) and (3)(4)(5)(6)(7)(8) are non-functional, vertex 6 -isolated.
Degradation of the FRC network through the edge removal is depicted in Table 3 (for constant boundary conditions) and in Table 4 (for gradient boundary conditions). The damage percentage in Tables 3,4  In the tables only non-isolated nodes are counted, that's why the number of vertices, inputs and outputs decreases. We also calculated the summary flow out of the LN (equal to the flow into the LN) and total flow through all system channels. In the tables we show a fraction of flow through the damaged system to the flow through initial undamaged system. The results depicted in Tables 3,4 reveal the following difference in the output of the flow for two types of pressure boundary conditions. In case of constant boundary condition total flow and outflow are higher, but the degradation rate is high as well. For gradient boundary conditions the flow continues even if 90% of system is destroyed. Pressure gradient on two neighbour boundary nodes allows a small flow between them, and it is not defined whether the boundary vertex is input or output. When 60% of the edges are removed the upper and lower halves of system loose their connectivity, but some small connected components remain. That's why there are less non-functional edges and more inputs and outputs in the «gradient» case than in «constant» case. Boundary points are highlighted in red. Figures 12,13 demonstrate the difference of pressure-coordinate dependance for graphs with consequently removed edges (pressure in isolated nodes is zero).

Graph Measures
From graph theory there are several observations and quantitative metrics that can be made of an FRC network that imply aspects of its function. Firstly, hierarchy, which can be defined as the imbalance between the number of nodes of a low degree and the number of nodes of a high degree. A network can have more lower degree nodes than higher degree nodes by adopting "hub and spoke" features where some higher degree "hub" nodes connect to many lower degree nodes. This is indicative of a scale-free network where the degree distribution is approximately that of a power law [41]. Currently, studies of FRC networks have found a Gaussian degree distribution indicating an absence of "hubs". A more formal definition can be used to explore this through two numbers; the mean local clustering coefficient, C, as defined by Watts-Strogatz and the mean shortest path length L. C i can be defined for each node, i, by the number of neighbours of i, N(i), which are also neighbours of each other, |N(i) ∩ N(N(i))| and the total number of neighbours, k(i) [42].
We can now define C as the average of C for all i, i.e. the average occupancy of connections between neighbours. If L(i, j) can be defined as the shortest possible path between two nodes i,j. Then L can be defined as, Where n is the number of nodes in the network. In lattice type networks the mean shortest path length is relatively long compared to random networks as is the clustering coefficient. Small-world networks are defined by their small mean shortest path length whilst still maintaining the high clustering coefficient of a lattice. In order to compare networks dimensionless forms of these numbers are defined.
Where the L ER and C ER forms of L and C are the quantities calculated on an equivalent Erdös-Rényi random networks. Erdös-Rényi equivalent random networks have the same number of edges on the same vertices, as the network, but with edges having an equal probability of existance between all vertices. From this the small-worldness quantity, σ can be defined as the ratio of the two [42].
For the network under consideration the quantities are as follows:Ĉ = 67.87,L = 0.817 ± 0.0605 and σ = 83.1 ± 6.19. The values suggest that this network has a similar mean shortest path length to a random network but a much greater degree of clustering giving a large small-worldness value indicative of a small-world type network. Analysis of actual FRC networks have shown average small-worldness values of 6.128 across 6 mice [14]. This is of particular relevance to FRC networks given their hypothesised role in maintaining chemokine gradients. A small-world network would be more likely to reduce gradients by spatially homogenising at "hubs" which will connect to a large number of "spokes". Maintenance of gradients between these spokes will be impeded by the "hub" compared to a lattice type network were in the absence of "hubs", a regular structure would exist between the "spokes". More complex gradients can exist across this regular structure than across a "hub", anode with a singular value. A greater understanding of the FRC network and the spatial variance of these quantities could allow more representative FRC networks to be generated.

Percolation Threshold
If each bond within a network has a probability, p, that it exists then the percolation threshold, p c , is the value for p below which a spanning tree cannot exist. That is the point at which the networks connectivity is sufficiently damaged that it is not possible to move across the network. At this point the largest spanning tree remaining is a fractal object which scales with the total size of the network according to a power law [43].
The percolation threshold for this network was found to be 0.32 with a standard deviation of 0.0825 which is close to the p c of a face centred cubic, 0.119. The network under consideration was subjected to a flow solver as described in Section 4.2 in order to study how fluxes change as p → p c . For each case the flux across the network was calculated, after which an edge was removed. This process was repeated until the flux was zero. 100 cases were run and the range and standard deviation for relative change in flux are shown below in Fig.14. As both viscous and diffusive effects have the same dependance on length the relationship shown in Fig.14 is true of both diffusion and fluid flux. The high p c suggests a high degree of topological robustness within the network implying a large difference in the importance of edges to maintenance of flow. The range indicates the minimum and maximum relative fluxes recorded at a given probability. Note that fluxes far below the standard deviation exist for all p > 0, which suggests that some edges are of considerable importance to flow. Verifying the existence of these properties in nature and correctly replicating them in artificial networks will be challenging.

Discussion
In this study we developed computational algorithms for modelling the geometry of two transport systems of the LN, i.e., the FRC network and the microvascular network. The first one provides the structure for lymph flow through the internal parts of the LN as well as cell migration. The blood vascular networks provide access for lymphocytes from the circulation to the LN parenchyma. As the remodelling of the networks takes place during infections, it is important to map the dynamic structure of the networks to their function in terms of cellular and information molecules (cytokines, antigens, etc.,) transport and distribution within the LN. The application of experimental imaging of internal LN structures is limited in humans. Therefore, computational modeling provides the tools for predicting the relationships between geometric and topological properties of the networks as well as on their performance under homeostatic and pathological conditions.
Capturing the 3-dimensional organization of LNs presents a challenge for existing imaging and analysis systems [21]. Recent studies provide basic information about geometrical and topological properties of the LN vascular systems. We used the data on organ-wide 3D-imaging of the microvascular network in murine LNs [21] to develop a voxel-based algorithm for 3D geometric modelling of the LN blood vascular network consistently with the data on blood vessels length distribution. A Cellular Potts Modelling framework was used to develop the FRC network model meeting the constraints on the volume and properties of FRCs. The blood microvascular-and the FRC network models were embedded into the SCS of the LN modeled as a sphere with the diameter about 200 µm. Taking together these two transport modules can be used to build up structurally complete computational models of the LNs as compared to those presented in [3,5,6] To our knowledge, this is the first study in which the lymph flow through the conduit network was studied under a range of physiological conditions. The FRC network can be severely destroyed during viral infection leading to immune deficiency [8]. It has been shown recently that an FRC network can tolerate a loss of approximately 50% of their FRCs without substantial impairment of immune cell recruitment, intranodal T cell migration, and dendritic cell-mediated activation of antiviral CD8 + T cells [14]. To evaluate the robustness of the conduit system in terms of the lymph percolation parameter, we have examined the lymph flow through the conduit network for an idealized geometry and under specified boundary conditions. The model based on a steady state Poiseuille equation predicts that the elimination of up to 60 − 90 % of edges is required to stop the lymph flux. This result suggests a high degree of the functional robustness of the network. We consider an idealised lymph node.
The conduits comprise of collagen fibers rather than being idealized empty cylinders [36]. We examined the impact of the presence of internal collagen fibers in the conduits on the resistance to fluid flow. Calculations, based on values in the literature of fibre size, conduit size and the typical number of fibres observed, seem to confirm that the collagen fibres inside the conduit are densely packed. This would suggest that the FRCs wrap themselves tightly around the collagen fibers. For a scenario of a conduit cylinder densely packed with collagen fibers, the computational flow study suggests that the diffusion would be the dominating process in mass transport as opposed to convective flow. These predictions warrant further experimental analysis. In addition, a more detailed analysis of our artificially generated FRC network showed that this network has a similar mean shortest path length to a random network but a much greater degree of clustering giving a large small-worldness value indicative of a small-world type network, which is qualitatively but not quantitatively matching to what was found experimentally [14]. In our study, we were interested in assessing the small-worldness property of the model network graph embedded into metric space, i.e. having realistic edge lengths distribution according to the real FRC network data. To this end we evaluated the shortest path length to be a real physical path length. Earlier, in [44] it was shown that the small-world parameter scales linearly with network size, for both model and real-world networks. The estimate of small-worldness parameter σ = 6.128 for murine FRC network data is based on the observed FRC network in T cell zone with about 170 nodes. In the model network for the entire lymph node the number of nodes is 3374, i.e., about 20-fold larger. The increase in the network size (the number of nodes) results in a proportional increase in the value small-world parameter, i.e. up to 83.1, consistently with the demonstration of the above cited study. These differences and their implications for flow and transport require further investigation. Finally, in an earlier study [45] it has been shown that fluid flow in turn regulates the reticular cell network, suggesting that models of the LN fluid homeostasis should consider elaborate feedbacks. LNs are characterized by complex multi-scale structures which represent critical issues in computational modelling of the LN physiology. Here we presented a consistent approach to integrating the blood microvascular-and the conduit network models within a confined space of the LN that will be used to study how the lymph-borne information is distributed to various parts of the lymph node under normal conditions and during an immune response. Further analysis should focus on studying the fluid-structure interactions for the vascular systems in conjunction with the cellular and chemical reactions taking place in live LNs. Finally, as it was stressed by W.E. Paul [46] "...It is to the quantitative prediction of the outcome of given perturbations in the immune system that we envisage our mathematical/modeling colleagues will apply themselves".