Computational Approach to 3 D Modeling of the Lymph Node Geometry

In this study we present a computational approach to the generation of the major geometric structures of an idealized murine lymph node (LN). In this generation, we consider the major compartments such as the subcapsular sinus, B cell follicles, trabecular and medullar sinuses, blood vessels and the T cell zone with a primary focus on the fibroblastic reticular cell (FRC) network. Confocal microscopy data of LN macroscopic structures and structural properties of the FRC network have been generated and utilized in the present model. The methodology sets a library of modules that can be used to assemble a solid geometric LN model and subsequently generate an adaptive mesh model capable of implementing transport phenomena. Overall, based on the use of high-resolution confocal microscopy and morphological analysis of cell 3D reconstructions, we have developed a computational model of the LN geometry, suitable for further investigation in studies of fluid transport and cell migration in this immunologically essential organ.


Introduction
The general functionality of the immune system is governed by a complex set of physical (transport), biochemical and biological processes occurring in space and time.Today, the need to embed immune processes into their spatio-temporal context is considered as a hallmark of systems immunology [1].The immune system consists of a network of interconnected lymphatic vessels, primary and secondary lymphoid organs (SLOs) of which lymph nodes (LNs) are the major sites for induction and maintenance of immunity.LNs have a highly elaborate structure-functional organization with various subcompartments orchestrating immune responses [2][3][4][5].The crosstalk between the hematopoietic and stromal cell compartments has a pivotal role in developing adaptive immunity in SLOs [6].One such stromal cell type are fibroblastic reticular cells (FRCs) which organize in a densely connected network and provide structural support for lymphocytes, secrete cytokines and chemokines, such as CCL19 and CCL21, which are crucial for proper lymphocyte migration [7].
Mathematical modeling of the spatial organization of the immune system and LNs in particular, is still in its infancy [8].Some studies have been performed in which the architecture of the LN is modeled as either 2D-or 3D regular orthogonal lattice [9][10][11][12][13][14] or a paradigmatic 3D view composed of a simplified set of basic structures [15].The broad availability of various imaging techniques has enabled the characterization of the LN structures with a high degree of spatial resolution.One of the first 3D reconstructions of the positioning and morphology of cortical sinuses, the T cell zone and high endothelial venules (HEVs) from confocal imaging data of an entire inguinal murine LN using Imaris (Bitplane) and MATLAB (MathWorks) software is presented in [16] for studying lymphocyte egress dynamics.The macro-structures such as B cell follicles and vascular network (HEVs) for popliteal murine LNs were reconstructed from OPT-based mesoscopic images analyzed with Imaris in [17] to study LN remodeling during a virus infection.Three major compartments of a mesenteric rat LN-sinus, paracortex and follicles-were 3D reconstructed in [18,19] for simulations of T cell transit through the LN.Recent 3D reconstruction and quantification of macro-structures such as B cell follicles, T cell zone and subcapsular sinus of an inguinal mouse LN from sections stained with immunofluorescence antibodies is based on original software (FOR3D) implemented in MATLAB [20].
The first model of the FRC network reconstruction was proposed in [21] using a static network of a 3D lattice that consists of thin, long rods with a random position and orientation.The in silico model of the FRC network defining the micro-architecture of the T cell zone was developed in [11], where the network edge length was based on imaging data from [22], and the motility of the T cell walk on the network was validated.To construct the FRC network, 3D lattices of nodes and edges were used with the FRCs randomly placed on a node in the lattice and the minimum number of fibers connecting the FRCs along the edges of the lattice set equal to two.In order to quantify the model parameters, two settings were considered: the "sparse" network structure based on the data from mice persistently infected with Lymphocytic Choriomeningitis Virus (LCMV) and the second "dense" (normal) setting with data provided from white pulp in the spleen of naive mice.In studies relevant for simulations of T cell movement [23], the reticular network was generated as a spatial random network without reference to topological data on FRC organization.In this study, a general mathematical model is described which is then parameterized with FRC edge length information from [22].The authors validate the random T cell motion on a network with this topology, generating quantitatively realistic migratory dynamics.
Computational generation of 3D morphology of biological tissues and organs using data from high-resolution microscopy has received increasing attention in recent years [24].Due to extremely variable histological appearance and resolution limitations of specific techniques, the generation of an entire LN from 2D scanned images seems to be infeasible yet.An alternative approach is to develop a generalized 3D solid model of the LN from a set of "elementary" macroscopic structures such as the subcapsular sinus, B cell follicles, trabecular and medullar sinuses, blood vessels, the T cell zone and FRC network.We sought to develop a LN model which included an anatomically realistic representation of the FRC network using high-resolution microscopy data, as FRCs have been shown to be vital for induction and maintenance of immune responses [6,7,21,22].We characterized experimentally the properties of the LN macro-structures and developed a computational method to reconstruct them in 3D.Finally, we generated a geometric LN model which is a crucial preliminary step necessary to establish an integrative LN model capable of further implementations of lymph flow, transport phenomena and cell migration.

Materials and Methods
In order to specifically target FRCs in murine LNs, a bacterial artificial chromosome (BAC)-transgenic mouse model was previously developed, which utilizes the Cre-loxP system under the control of the Ccl19 promoter [25], referred to as Ccl19-Cre x R26-eyfp reporter mice.The Ccl19-Cre transgene activity targets past and present CCL19 production in LN FRCs in vivo.
A representative immunofluorescence Z-stack image of the murine LN is shown in Figure 1a.The major visible structural units are the subcapsular sinus (SCS) stained by α-smooth muscle actin (SMA), efferent lymphatic vessel, B cell follicles (B220) and the FRC network (EYFP).For the geometric modeling, we quantified the relevant sizes as follows: the LN diameter ~500 µm, the arteries and veins diameter 15 µm.The B cell follicles can be viewed as asymmetric biconvex lenses having a diameter of 150 µm and height of ~75 µm.Lens refers to the shape of a 3D object obtained by rotating a 2D lens about its narrow axis of symmetry.
The most complex geometric structure in the LN is the FRC network as shown in Figure 1b.The area shown is a representative paracortex T cell zone which is occupied by FRCs, sufficiently large for quantitative morphological analysis.The T cell zone FRCs themselves occupy approximately 3%-4% of the total volume, estimated by 3D reconstruction of the EYFP+ FRC network from n = 7 mice using Imaris (Bitplane).The diameter of an FRC body is 5-7 µm and the average distance between centers of mass of FRCs is ~23 µm, determined by isolating single FRCs from 3D reconstructions in Imaris.Each FRC contains 4-5 protrusions (median) and the protrusions have a very variable length.The overall summary of the statistical distributions of the degree and the length of the edges for the FRC network is shown in Figure 2.

Results and Discussion
The visualization of the solid objects in the computational LN model was performed using 3D computer graphics software Autodesk 3ds Max.The solid elements were constructed using original algorithms written in C++ language implemented in Microsoft Visual Studio 2012.

B Cell Follicles, Trabecular Sinuses, Blood Vessels and Medulla
B cell follicles were reconstructed using LN confocal data.In the homeostatic case, which we consider, the B cell follicles are rather stable in forming spheroid-type shapes.The basic starting structure of the B cell follicles is defined a priori using original data quantitatively consistent with an independent study [20].By marking several points along the B cell follicle, we created a piecewise-linear approximation of the B cell follicle border as shown in Figure 3a.Spline interpolation was used to generate a smoother model approximation [26].The next step was to set up the number of triangles on the top of the solid object (Figure 3b), which is a variable parameter that allows one to control the resolution of the B cell follicle discretization.The overall algorithm can be summarized as a pseudo code using the following stages: (1).
Initial piecewise-linear approximation of the B cell follicle border is created by marking several points on the confocal image.(2).
Generating an updated set of the border points.(4).
Setting up the number of triangles on the surface of the object to achieve a required mesh resolution to be used for computational discretization of the B cell follicle.
Trabeculae and the medulla area were modeled using a piecewise-linear approximation as well.Trabeculae are a mix of connective tissue and small amounts of plain muscle fibers, penetrating into the LN and they were constructed separately from the SCS inner membrane.A piecewise linear trajectory starting from the membrane and consisting of vectors orientated in different angles was used to create a jogged line.In the next step this line was converted into a 3D object and random numbers were added to the coordinates of all points on the surface to obtain the final object as shown in Figure 3c.In general, the term "medulla" denotes the entire macrophage-rich, loose region around the efferent lymphatics.In the above computational construction we restrict the notion of "medulla" to just the efferent lymphatic vessels.Consideration of the extended boundaries of the "medulla" domain can be incorporated within the solid modeling approach once the quantitative data on the material properties are available.The blood vessels and high endothelial venules (HEVs) in the LN were approximated using the jogged line approach as described above and then triangulated (Figure 4).In general, the appearance and localization of the trabeculae, blood vessels and cortical and medullar sinuses is highly variable between different LNs [16,20].The constructions of trabeculae, blood vessels and medullar sinuses are based on our original 3D confocal data including those from [17].Using the developed code, the geometric model can be further tuned to take proper account of specific data on shapes and localization of trabeculae, blood vessels and medullar sinuses.The outer and inner surfaces of the SCS were approximated as a triangulated sphere with a noise factor in the mesh coordinates (Figure 5a).The noise was used to make the appearance of the SCS closer to a realistic LN.Indeed, the boundary of the SCS, as shown in imaging-based reconstructions [20], has a different appearance than an ideal sphere.It is rather an irregular surface with mathematically unknown yet pattern properties.As the structural characterization of the SCS shape remains poorly analyzed, we follow the principle of the Central Limit Theorem, indicating that the most reasonable choice for the distribution of a random variable is Gaussian, in the absence of any other information [27].Therefore, we used normally distributed random variables with zero mean and the standard deviation of 1.6 μm to perturb the mesh coordinates of an ideal sphere, thus representing more realistic SCS geometry.In addition, the perturbation was checked for not being too large for the nodes with small coordinates.The integrated view of the reference LN with macroscopic solid objects constructed in this section is shown in Figure 5b.

Generation of the FRC Network
The process includes the following steps described below: (1).Generate the FRC network tree according to the experimental statistical data.
(2).Build the voxel approximation according to the FRC network tree.
The computational generation of the FRC network includes the following stages.The first step is creation of the FRC network tree.The first mesh node is placed in the origin, i.e., at the point with coordinates (0, 0, 0) and in the next step, vectors (Xi, i = 1,…,n) are generated according to the parameters specified in Table 1 for the lengths of edges between FRCs and number of edges (degree) for single FRCs.Each vector group is checked whether the sum of the vectors is close to zero ("near zero value" check): for vectors (X1,…, Xn) it computes the mean value Xs = [sum(Xi)/n].The group is generated and checked again if Xs is greater than defined parameter [possErr].The final point [Pf] of each generated vector can become a next FRC node if the distance between [Pf] each of the existing nodes is more than [minDist] parameter.Otherwise, a new node will be randomly generated.This process continues until one of the edges reaches some specified object (e.g., the surface of the B cell follicle).The second step is a voxel-based approximation of the FRC tree created as described above.Each voxel of the tree is a cube with specified edge length (voxLength).The algorithm computes the centers of the two adjacent voxels as follows: where (x1, y1, z1)−>(x2, y2, z2) are the coordinates of some subvector of tree vectors.Each tree vector is represented like an array of subvectors with length less than the voxel size.
The voxel-based approximation of the FRC network is shown in Figure 6a.As the FRC network is localized to the LN paracortex, it is necessary to add spatial restrictions while constructing the geometric model.The algorithm selects two vectors of each polygon, normalizes them according to the size of the voxel and generates a set of linear combinations that do not exceed the boundaries of the designated domains.New voxels are then generated according to the calculated set of points.After the third step, the overall set of voxels consists of three voxel groups: (1) FRC network voxels reachable (connected to) from the first node at the origin (0, 0, 0); (2) FRC network voxels non-accessible from the first node and (3) a set of voxels that touch the domain border.The algorithm identifies and creates all the voxels of the FRC network that can be reached from the first node (0, 0, 0).Some voxel pairs could touch via edges or vertices and they are considered as incorrect segments in the set.The algorithm fixes them by adding new voxels around the problematic contact zone.The extension of the voxel set continues iteratively until no incorrect segments are found.
The fifth step of the algorithm generates a pre-final solid model of the FRC network.It converts voxel coordinates to real coordinates in μm and each voxel is represented as a set of 12 triangles.Importantly, only those faces that do not touch the faces of neighboring voxels are selected for the model.Each vertex is also checked for its uniqueness.The problem of enumerating the voxels without duplications is a computationally demanding task of the algorithm.
Finally, the algorithm performs the smoothing of the voxel-based structure to achieve biologically appealing shapes of the FRCs.For each vertex [Vc] the algorithm calculates the midpoint [Vm], which depends on data from the pre-final solid model and vector Vn = [Vc, Vm].The vertex is then assigned a new coordinate Vc * = Vc + [alpha]*Vn, where [alpha] is a smoothing parameter in the range (0, 1).The result of the smoothing of the local FRC network is shown in Figure 6b, morphologically consistent with the appearance of the FRC network in vivo (Figure 1b).Further modifications of the FRC network can be implemented once additional topological data become available.

Generation of the Whole LN
The final 3D geometric model of the LN is shown in Figure 7a.It consists of a set of embedded objects and has the overall complexity specified by the following parameters: As a result of retriangulation of all objects, the LN geometry was reproduced with 621,376 triangles and the FRC network required 1,759,180 triangle polygons.For Intel Core i7-2670QM CPU @ 2.20 GHz the computation time of the LN model generation, namely the assembly of all macroscopic elements, is approximately 1 h.It is clear that the overall complexity of the FRC network is much more demanding as compared to the other macroscopic LN structures.The consistency of the FRC network with the biological data can be appreciated from the 20 µm slice of the LN model (Figure 7b) and the following statistics: the average number of degrees per node (links per FRC) is 4.2, the length of edges ranges from 0.24 to 44.73 µm, with a mean value of 17.5 µm, the subcapsular sphere with diameter of 500 µm contains ~10 4 nodes (FRCs).The source code and the LN model data file in the below specified formats are available upon request.We analyzed the statistical properties of the computationally constructed FRC network in order to validate the model compared to experimental data.The distribution of the number of edges per node and the length of the edges is shown in Figure 8.
The solid LN model can be further processed to generate an adaptive mesh model of the LN suitable for finite-element discretization of the transport equations.The tetrahedral mesh approximation of the B cell follicles, the trabecular and medullar sinuses of the LN model is shown in Figure 9, generated using the open-source adaptive mesh generator Ani3D (http://sourceforge.net/projects/ani3d/).A number of formats are supported, including Ansys, Acis, STL, Gmsh, OpenFOAM, OBJ.

Conclusions
This is a first study in which a 3D solid model of the LN geometry has been systematically developed in a biologically relevant manner.A computational algorithm has been developed to generate the macroscopic structures of a murine LN based on high-resolution confocal microscopy data.The technology provides a bridge between existing models of immune processes based on a simple 2D-or 3D lattice representation of the LN morphology and anatomically realistic LN models.A library of solid objects modeling the geometry of SCS, B cell follicles, trabecular and medullar sinuses, blood vessels, the T cell zone and FRC network has been created.The models can be exported in any format, e.g., OBJ, STL, making it suitable for further use in computational studies of lymph and cytokine transport/distribution and chemokine-guided cell migration in the LN.
Overall, a reference geometric model of the 3D LN structure consistent in its resolution with the generalized imaging data-driven view of LN organization has been developed.The model can be further adapted for individual LNs once specific structural data are available.In the present approach, the basic starting structure of the B cell follicles has to be defined a priori; it cannot emerge from another property.The whole LN structure changes in the course of an immune response, but the underlying processes are beyond the scope of our study.The generated meshed model of the LN in conjunction with the computational algorithm implementing the 3D model represents an important step towards studying transport phenomena in SLOs and the dynamics of LN remodeling during viral infections is a further milestone we aim to computationally model in the future.
The 3D geometry of LNs is highly complex and a mathematical model of LN architecture is a necessary pre-requisite to computational fluid dynamics studies.The 3D solid model of the LN can be used via adaptive 3D meshing and finite-element discretization of the governing equations in computational studies of a broad range of transport phenomena in LNs, including the prediction of cytokine fields, chemokine gradients and pharmacokinetics of drugs [28].
It has been recognized that space is the current frontier in immunology.Capturing the structure-function relationship in a computationally efficient manner is the key to successful systems immunology approaches [29].It requires anatomically-based mathematical models integrating interaction processes across multiple scales.This sets an agenda for further interdisciplinary research on SLOs, which should address the following aspects: the material properties, blood and lymph flow patterns, fluid-tissue interactions, and last but not least -the fine architecture of the lymphoid organs.The results of our study clearly indicate that the modern technologies available in computational mathematics enable a systematic advance of systems immunology to a genuinely integrative approach, which can be further utilized in studies of regulation of the immune system in health and disease.

Figure 1 .
Figure 1.(a) Lymph node (LN) architecture with major structural units: subcapsular sinus (red), efferent lymphatic vessel (red), B cell follicles (cyan) and the FRC network (green).White rectangle indicates representative T cell zone; (b) High resolution Z-stack of the topology of the T cell zone fibroblastic reticular cell (FRC) network.Podoplanin (PDPN) is also used as a marker for FRCs.Scale bars represent 200 and 30 μm, with Z-depth 38.6 and 23.6 μm, respectively.

Figure 2 .
Figure 2. Fundamental characteristics of the FRC network: the nodes represent the FRC body and edges are physical links between them.The edges per node-and the edge length distributions for the T cell zone FRC network.Data obtained from n = 7 mice from two independent experiments.

Figure 3 .
Figure 3. (a) 2D spline-based approximation of the B cell follicle; (b) 3D approximation of the B cell follicle; (c) 3D solid model of the trabecular sinus; (d) 3D solid model of the medulla.

Figure 4 .
Figure 4. Basic solid models for (a) blood vessels and (b) vessels branching in the high endothelial venules (HEV) domain.

Figure 6 .
Figure 6.Voxel-based generation of the FRC network.(a) Initial local structure; (b) Smoothed solid model of the local structure.

Figure 7 .Figure 8 .
Figure 7. 3D geometric model of (a) the whole LN and (b) the 20 µm slice showing the FRC network, B cell follicles and the medulla.

Figure 9 .
Figure 9. Adaptive mesh for B cell follicles, the trabecular and medullar sinuses of the geometrical LN model generated using tetrahedral spatial discretization software Ani3D.

Table 1 .
The parameters of the algorithm for generating the fibroblastic reticular cell (FRC) network.