Bio-physically plausible visualization of highly scattering fluorescent neocortical models for in silico experimentation

Background We present a visualization pipeline capable of accurate rendering of highly scattering fluorescent neocortical neuronal models. The pipeline is mainly developed to serve the computational neurobiology community. It allows the scientists to visualize the results of their virtual experiments that are performed in computer simulations, or in silico. The impact of the presented pipeline opens novel avenues for assisting the neuroscientists to build biologically accurate models of the brain. These models result from computer simulations of physical experiments that use fluorescence imaging to understand the structural and functional aspects of the brain. Due to the limited capabilities of the current visualization workflows to handle fluorescent volumetric datasets, we propose a physically-based optical model that can accurately simulate light interaction with fluorescent-tagged scattering media based on the basic principles of geometric optics and Monte Carlo path tracing. We also develop an automated and efficient framework for generating dense fluorescent tissue blocks from a neocortical column model that is composed of approximately 31000 neurons. Results Our pipeline is used to visualize a virtual fluorescent tissue block of 50 μm3 that is reconstructed from the somatosensory cortex of juvenile rat. The fluorescence optical model is qualitatively analyzed and validated against experimental emission spectra of different fluorescent dyes from the Alexa Fluor family. Conclusion We discussed a scientific visualization pipeline for creating images of synthetic neocortical neuronal models that are tagged virtually with fluorescent labels on a physically-plausible basis. The pipeline is applied to analyze and validate simulation data generated from neuroscientific in silico experiments.


Background
Scientific visualization is a key component in neurobiology. It helps neurobiologists to explore and convey different levels of interpretations of complex sets of neuroscientific data. Recent advances in computational sciences and hardware technologies allowed some biological experiments to move from the wet laboratory to computer simulations, to in silico [1][2][3] experiments.
This paradigm shift is expected to accelerate and consolidate the research discovery and also to enable novel capabilities in the near future. It will reduce the dramatic costs of clinical trials and complement the traditional in vivo and in vitro methods [4,5]. Nevertheless, this approach requires developing rigorous mathematical models of the biological experiments and their surrounding physical conditions and then plugging them in high performance computer simulation applications. These applications are designed to exploit the growing computing power of stateof-the-art supercomputers to simulate and analyze complex biological processes at different scales of resolution [6].
This emerging trend opens novel avenues for multiscale computational modeling of the brain tissue, and in turn a better understanding of how the brain works.
In this context, visualization is not merely exploited for providing visual analysis of the data; it is a significant tool for evaluating and validating the results of in silico experiments. This visual feedback closes the loop and affords the neuroscientists an effective environment to tune and enhance their models and also to improve the accuracy of the simulations in an iterative manner.

Motivation
The current neuroscientific visualization tools have been improved considerably during the last years to visualize simulation data. A clear example is given by Hernando et al. to interactively visualize the simulation of the cortical activity of large scale neuronal microcircuits [7]. Nevertheless, the toolset is still inadequate for visualizing and validating the data generated from various in silico experiments such as voltage sensitive dye imaging (VSDI) [8], Calcium imaging [9] and also optogenetic stimulation experiments. For example, visualizing the data arising from simulating an optogenetic procedure entails incorporating plausible optical models into the visualization pipeline to account for light interaction with highly scattering turbid media [10]. Accurate visualization of the responses from simulated imaging experiments requires a sophisticated bio-physically-based optical model that incorporates fluorescence in the rendering integral and can account for the actual optical properties of the biological tissue. Such pipeline is still largely unfulfilled and will require an extensible spectral visualization system that can model and simulate light interaction with highly scattering fluorescent volumetric data resembling the fluorescent structures in real tissue.
We address these shortcomings and present an advanced visualization pipeline that can accurately render highly scattering fluorescent volumetric datasets. This pipeline is mainly applied to a fluorescent brain model that represents a digital reconstruction of the microcircuitry of somatosensory cortex of rats [9] to validate its structural and functional aspects. For instance, it is currently used to perform in silico VSDI experiments for validating the cortical activity of the reconstructed model against in vivo imaging experiments [8]. Moreover, it can be useful for other fields such as computational microscopy, where a physically-plausible simulation of microscopic fluorescent images is required for analysis purposes [11][12][13].
Our pipeline is composed of two software workflows. The first one is a generic physically-based visualization engine for rendering highly scattering heterogeneous fluorescent volumes. The other workflow is developed in particular to efficiently extract a fluorescent tissue block volume from the neocortical column micro-circuit presented by Markram et al. [9].

Design and implementation of an extensible pipeline
for visualizing fluorescent-tagged scattering volumetric datasets. 2. Rigorous physically-based optical model to simulate the light interaction with fluorescent participating media, taking into account their spectroscopic and optical properties. 3. Qualitative validation and analysis of the developed optical model by correlating the spectral power distributions (SPDs) (or responses) of the generated images with respect to experimental emission spectra of different fluorescent dyes. 4. Design and implementation of an automated parallel workflow for generating an extracted fluorescent tissue block from the neocortical column model. 5. Visualization of fluorescent neuronal models tagged with multiple fluorescent solutions having different optical properties and evaluating the results collaboratively with neurobiologists.

Related work
Neurobiology scientists are familiar with generic visualization packages such as Paraview [14], Voreen [15] and ImageVis3D [16]. They use them frequently to visualize and analyse data acquired from sensing devices, for example imaging scanners and microscopes. In some cases, these software packages can be employed for visualizing certain structural aspects of the data arising from in silico experiments and modeling procedures, for example, to validate the morphological distribution of the neurons in the neocortical column model [9]. Other frameworks have been customized to fulfill specific demands required by the scientists such as Voxx [17] and VAA3D [18]. The design goals of the previous frameworks have been focused on scalability and interactivity. Consequently, they traded the performance with oversimplified optical models that remain very limited to visualize fluorescent data or even to enhance the photorealism of the generated image [19,20]. Photorealistic visualization of neuroscientific data with advanced illumination models was addressed in two studies. The first one is presented by Banks et al. [21]. They integrated global illumination into their visual data analysis pipeline for displaying the fiber tracts of the brain. Their study was intended to improve the data interpretation in the presence of complex jungle of fibers surrounding brain tumors. The other study presented Exposure render, an interactive GPU-based framework that coupled Monte Carlo ray tracing with physically-based light transport models to generate highly realistic renderings of volumetric data [22]. This framework is capable of visualizing in silico optogenetic experiments, but it cannot be employed to visualize fluorescent data.
Visualizing fluorescent volumetric data was firstly presented in FluVR [23], a commercial application that used a simple deterministic physically-based model called the simulated fluorescence process (SFP) to combine elastic and inelastic rendering. Although it was capable of handling multiple fluorescent dyes in the volume, FluVR was limited in several regards. The SFP assumed that the emission occurs only at the maximum emission wavelength and ignored the rest of the emission spectrum. This optical model did not account for the spectral characteristic of the dyes and ignored multiple scattering.
Physically-plausible visualization of fluorescent participating media has been investigated in few computer graphics research studies. These studies were exclusive to specific applications and their implementations were not developed in the form of an integrated framework that could be utilized for other purposes. In summary, these studies have developed extensions to integrate the fluorescence phenomena into the rendering equation [24][25][26][27][28], but they were limited to certain extent. Glassner [24] presented the first formalism of the full rendering equation to simulate the fluorescence effect. However, his model ignored the distinct properties of the fluorescent dyes. Cerezo [27,28] and Gutierrez [25,26] have extended Glassner's model to account for these missing parameters. Nevertheless, their models used biased rendering methods (discrete ordinates and curved photon mapping) to render the fluorescent pigments of the ocean. Moreover, they ignored the actual spectral properties of the dyes and used oversimplified profiles for the excitation and the emission spectra. Abdellah et al. presented a physicallybased framework for simulating imaging experiments with light sheet fluorescence microscopy. The optical model developed in this study presented further extension to the previous fluorescence models taking into account the intrinsic characteristics of fluorescent dyes [29,30]. They also validated their model against realistic emission spectra of multiple fluorescent dyes. This model was only capable of visualizing tissue models with negligible scattering properties to simulate the imaging of clarified brain tissue [31,32], but it failed to handle volumetric tissue models with highly scattering content. Our optical model presented in the following section is introduced to fill this gap.

Optical models
Based on ray tracing and the basic principles of geometric optics, advanced optical models of volume rendering ideally solve the radiative transfer equation (RTE) to simulate the light transport in a continuum and generate a physically-plausible synthetic image [33][34][35]. The general formulation of the light transport presented by Veach [36] is extended by Pauly et al. [37] to handle scattering media. Nevertheless, this formulation has never been investigated for considering the fluorescence effects. In the following part, we begin with this extension to derive the path integral formulation of our fluorescence optical model. Table 1 summarizes all the relevant terms and symbols that appear later in the text. We also recommend the reader to refer to [38] for further explanation of some of the terms in the rendering integrals.

Path integral formulation in fluorescent volumes
Assuming a path consisting of three points x 0 x 1 x 2 , where the light source and the camera film are located at points x 0 and x 2 respectively, and x 1 is sampled to be a random interaction point in the volume where the light scattering occurs (Fig. 1), the radiance arriving to the camera following a scattering event at x 1 can be computed with the monochromatic light transport formula described in Eq. (1), where ω = x 0 ← x 1 is the incoming direction, L ve and L s are the radiance due to self emission and scattering respectively. The self-emission term is usually ignored unless the volume itself is emitting due to chemical or thermal processes, which is out of the scope of the presented model. In this case, the total radiance recorded by the camera due to light scattering L s in the volume is evaluated with the integral in Eq. (2), where σ s and f p are the scattering coefficient and the phase function of the volume respectively and L i is the incoming radiance towards the point x 1 from any direction ω .
For convenience [39], Eq. (2) can be re-written in the form of Eq. (3) as an integral over surfaces dA and volumes dV instead of directions dω on the sphere 4π to yield what is called the three-point form of the light transport equation, where F s , G,V , V , τ and L e are the scattering function, geometric term, visibility term, binary visibility function, transmittance and the emitted radiance from the light source at x 2 respectively.
If the light scatters at n−1 interaction sites before reaching the camera at x 0 , where x n is a sampled point on the light source, the path integral equation becomes where L e is the emitted radiance from the light source at the sampled point on its surface x n to the first interaction point in the volume x n−1 .
In principle, Eq. (8) can be used to render highly scattering volumetric models assuming monochromatic wavelengths, i.e. there is no transfer of energy from one wavelength to another. We have extended this equation by introducing a term called the path binary fluorescent visibility V f i that indicates whether a path has encountered a fluorescence emission or not. Adding this term to Eq. (8) and integrating over all excitation wavelengths λ x to evaluate the radiance at specific emission wavelength λ m , the rendering equation becomes

Monte Carlo estimator
The path integral formulation of our fluorescence model, Eq. (9), evaluates the radiance arriving to the camera at point x 0 from direction ω at certain emission wavelength Fig. 1 Light transport in a highly scattering volumetric extent. a The volume prior to illumination by the light source. b Single scattering interaction: the light ray is scattered once between the light source and the camera on a single path x 0 x 1 x 2 . c Multiple scattering: the light ray bounces multiple times between several interaction events before reaching the camera on a single path x 0 x 1 x 2 . . . x n−1 x n . d The radiative transport equation evaluates the light propagating from the light source to the camera on multiple paths x 1 , x 2 , . . . , x n . The rays are shot from the camera towards the light source to sample the scattering events λ m after multiple scattering events in a highly scattering fluorescent volume. In a stochastic path tracer, this integral can be approximated with the Monte Carlo estimator expressed by Eq. (10), where p(.) is the probability density function (PDF) for sampling a point x n on the surface of the light source, an excitation wavelength λ x from the emission spectrum of illuminating light, a scattering event with a direction ω j and a distance t j . The path binary fluorescence visibility term V f i accounts for the spectral optical properties of the volume, the intrinsic spectroscopic properties of the fluorescent dye including its excitation and emission spectra, molar absorptivity and quantum yield, and also the concentration of the fluorescent solvent in a given solution. where Monte Carlo path tracing is used to determine the interaction sites, or events, within the volume extent. The fluorescent events -represented by the green points in Fig. 2 -are stochastically identified according to the ratio between the fluorescence absorption coefficient μ f a and the total absorption coefficient μ a of the volume at emission wavelength λ m . There are eight possible combinations that might occur during the path sampling. According to the type of the sampled event, some of these cases are plausible and the other are not possible as explained in Fig. 3. The SPD of the fluorescence absorption coefficient μ f a (λ) is expressed in terms of the excitation (or absorption) spectrum of the fluorophore f x (λ), the concentration of the dye in the solution C, and its molar absorptivity at the maximum excitation wavelength . The spectral radiance is computed by tracing a ray through the volume at certain wavelength between 300 and 800 nm with 1 nm increments. The estimated pixel value is updated only if the constructed path is valid and a fluorescence emission occurs. A valid contributing path, such as 2 and 4 in Fig. 2, consists of a series of elastic scattering events and a single inelastic one that involves changing the wavelength from λ x to λ m . In this case, the light source is sampled and the radiance emitted towards the fluorescence emission event is attenuated according to λ x . Otherwise, the fluorescence visibility V f term is set to zero and the path is terminated. The paths are sampled with woodcock tracking, which is known to be an unbiased method [40,41].
The probability of fluorescence emission p f is expressed by two terms: the photon absorption probability p x and the photon emission probability p m [42] Therefore, the fluorescence emission probabilistically occurs in terms of the exact spectral characteristics of the fluorescent dye including its excitation f x (λ) and emission

Virtual fluorescent tissue volume generation
The digital model of the neocortical column is organized in a circuit, which can be seen as a database containing a set of neurons having diverse morphological and electrical characteristics. These neurons are statistically positioned and oriented within the 3D extent of the column [9,43]. The virtual fluorescent tissue block is reconstructed from the neocortical column circuit for in silico experiments in five basic steps (Fig. 4): 1. Identifying a list of neurons that will be contained in the resulting tissue block. This list can be selected based on common morphological or electrical properties to address specific kind of in silico experiment. 2. Creating a watertight surface mesh model for the block from the morphological descriptions of the neurons in the circuit. If a given morphology is broken, the neuron identifier is reported to fix the morphology. The neuronal morphologies are converted into watertight surface meshes using an extended version of the workflow presented by Lassare et al. [44]. The individual meshes generated for every neuron are loaded into Blender [45] and the final mesh block is extracted based on the extent of the requested block.
3. Converting the mesh model to a volumetric one using solid voxelization [46]. This operation is handled with a fast in-house GPU-based voxelization software that uses conservative rasterization [47]. If the input mesh is not watertight, the neuron identifier of the mesh is reported to be fixed. 4. Annotating the volumetric tissue block with the optical properties of the rat brain at the specified region. The optical properties are retrieved from a 3D atlas that was compiled in a recent study by Azimipour et al. [48]. 5. Labeling the block with fluorescent dyes to simulate their injection into the intracellular space of the different neurons contained in the generated block. The intrinsic spectroscopic characteristics of the selected dyes are obtained from an online database available at [49].
In some cases, the experiments are limited to investigate the responses of individual neurons, pair of neurons or a small set of neurons. The generation of a fluorescent tissue block for such experiments is relatively easy as described in the aforementioned process. In contrast, other experiments require extracting a large tissue block that might assemble hundreds or thousands of neurons. The spatial extent of this block does not necessarily enclose the bounding volumes of all the neurons that are located into it because the positions of the neurons are identified based on their cell bodies (or somata). Extracting a tissue block from a large cluster of neurons following the previous approach on a single computing node is inefficient and in some cases is impractical. To resolve this issue, we have developed a parallel workflow that can efficiently generate high density tissue blocks. This workflow runs on high-end visualization clusters that consist of several computing nodes connected together via high bandwidth networking infrastructure. This workflow, shown in Fig. (5), parallelizes the mesh generation and clipping operations exploiting all the available nodes in the cluster.

Pipeline implementation
Implementing our optical model requires a physicallybased spectral rendering framework that can model the light rays by spectral distributions as an alternative to the tri-stimulus representation. The physically-based rendering toolkit (PBRT) [50] has been chosen amongst other systems like Mitsuba [51] or LuxRender [52] due to the existence of an accompanying reference [38] that  5 An illustration of the mesh block extraction from the selected targets in the cortical column. a The spatial extent of the block is identified by a bounding box that is given in the input configuration. b The meshes are generated from the corresponding morphologies with an extended version of the meshing pipeline presented by Lasserre et al. [44]. c The resulting wavefront object meshes are loaded in Blender [45] and clipped on a per-mesh basis. d All the clipped meshes are loaded in Blender and grouped together with a union boolean operation to generate the final mesh block documents the software architecture of the framework. Though, it only supports CPU-based rendering, which will limit the rendering performance for high resolution images with sufficient sampling densities.
We have implemented our estimator in Eq. (10) in a volumetric integrator class that can be selected in the configuration file given to run the rendering framework. We have also extended the volumetric grid class to support loading annotated fluorescent volumes to allow tagging the same model with multiple fluorescent dyes. The automated block extraction pipeline is configurable to generate PBRT scene description files and render them directly after the creation of the fluorescent tissue block volume.

Results, validation and discussion
The results of our visualization pipeline are demonstrated on a 50 μm 3 tissue block extracted from the center of the neocortical column model (Fig. 4). A surface rendering image of the surface mesh of this block (prior to virtual fluorescent injection) is illustrated in Fig. 6.
From this extracted mesh block, we have created two experimental sets of fluorescent-annotated volume blocks. The first one is tagged with the same type of fluorescent dye dissolved in several solutions having different extinction coefficients. The goal of this set is to experiment the responses of the same fluorescence parameters in the presence of relatively low, medium and high scattering volumes. The other set is labeled with various fluorescent dyes that have different spectral responses at fixed concentrations. This set is designed to validate and measure the performance of our extended optical model that can simulate the light interaction with fluorescent volumes. The two sets were labelled with multiple dyes from the Alexa Fluor family, Alexa Fluor 350, 488, 568 and 633. This family is selected in our experiments due to its importance in fluorescence microscopy and cell biology in general [53]. Table 2 summarizes some of the spectroscopic properties of the four dyes including their maximum excitation and emission wavelengths (nm), molecular weight (kDa), quantum yield, and molar absorptivity (cm −1 M −1 ).
The first set is labelled with three Alexa Fluor 488 solutions that are characterized with extinction coefficients that are 10, 100 and 1000 times greater than that of pure  water [54]. To maximize the emission, the illuminating light source is set to emit at the maximum excitation wavelength of Alexa Fluor 488 at 495 nm. Figure 7 shows the results of rendering the three tissue volume blocks under the same illumination conditions. The tissue blocks in the second set are tagged with Alexa Fluor 350, 488, 568 and 633 solutions at the same concentration (0.4 mol/l). The same illumination conditions defined in the first experiment are used to excite the volumes in this case where the light source emits at the maximum excitation response of each respective dye (refer to Table 2). Figure 8 illustrates the images rendered for the four tissue volume blocks used in this experimental set.

Fluorescence optical model validation
The experimental measurements of the excitation and emission spectra of fluorescent dyes are normally recorded for highly diluted and low scattering solutions using Beer-Lambert law and the fluorescence brightness equation [55]. However, the normalized spectral distributions of the emission spectra recorded from highly scattering solutions should have similar profiles to the experimental emission spectra of the fluorescent dyes [56]. In this context, we validated our fluorescence optical model relying on two basic tests. The first one measures the SPD of the generated images from our visualization pipeline and then compares their normalized profile with the distribution of the intrinsic emission spectra of each dye. Note that the SPDs of each image are recorded before their conversion to RGB colors for each pixel in the image.
The four tissue volume blocks in the second experimental set are used to validate our optical model. The normalized spectral responses (or SPDs) from the four images shown in Fig. 8 are compared to the emission profiles of the four dyes. The results of this validation test are shown in Fig. 9.
The second validation test measures the performance of the model when the volume is illuminated with different wavelengths. Depending on the excitation spectrum of the dye and the selected wavelength to illuminate the solution, the scale of the emission spectrum is proportional to the amplitude of the excitation spectrum at the excitation wavelength. The maximum emission profile is reached when the maximum excitation wavelength is used [55,57]. In this test, all the tissue volume blocks are illuminated at several wavelengths (300, 346, 495, 532, 555, 578, 632 and 700 nm) and the responses are recorded and relatively compared. The results of this test are illustrated in Fig. 10.

Pipeline evaluation
The rendering results were evaluated collaboratively with a group of different experts in neurobiology and in silico neuroscience. They all agree that the renderings are similar to what they visualize under the microscope. They were also excited to see how the responses are changing when the optical and spectroscopic properties of the dyes are varied. This would allow them characterizing the responses of the neurons in various regions of the brain that have different optical properties. The scientists  working in the brain simulation team have expressed their interest in applying our pipeline to their data to validate their in silico VSDI experiments against realistic data recorded by the fluorescence microscope. Other scientists have requested further extensions of the pipeline to visualize neuroglial cells.

Rendering performance
In general, Monte Carlo rendering is known to be time consuming. The rendering performance of Monte Carlobased algorithms depends on multiple factors including the pixel sampling density, number of light samples, optical properties of the volume and the image resolution as well. If the sampling rates are relatively low, the rendered image will be full of noise. Therefore, high sampling is mandatory to have an image with a converging solution. Our results have been rendered with pixel sampling of 512 × 512 samples per pixel. Moreover, high spectral sampling is also required to obtain accurate emission spectra that can reflect those measured in real spectroscopic experiments. We have used a spectral sampling of 1 nm. The rendering time of the images demonstrated Fig. 9 Normalized emission SPDs measured from the images illustrated in Fig. 8. The spectral responses of the emission recorded from each tissue block is qualitatively compared with the actual emission spectra of the four Alexa Fluor dyes used to tag the tissue block. The SPDs are obtained at the maximum excitation wavelengths of each respective dye (346, 495, 578 and 632 nm) and 1024 spectral samples per pixel

Conclusion and future work
The current visualization systems are limited to meet the immense challenges of in silico neuroscience, where biological experimentation are performed in computer simulations. A wide range of those experimental observations rely on fluorescence imaging to reveal several structural and functional aspects of the brain. Reproducing the same experimental procedures in silico is subject to the existence of visualization engines that can handle fluorescent models. We presented a visualization pipeline to address these challenges. The pipeline is composed of a generic volume rendering system capable of visualizing highly scattering fluorescent volumetric datasets. This system is applied to visualize virtually-tagged fluorescent tissue blocks that are extracted from a unifying model of the neocortical microcircuitry reconstructed from rats. The pipeline is primarily developed to assist the neuroscientists exploring and analysing their in silico experimentations that incorporate those fluorescent blocks to present a visual feedback that allows them fine tuning their experimental parameters and improving the model in an iterative manner (Fig. 11).
A rigorous bio-physically-based optical model is developed to account for light interaction with highly scattering fluorescent media. This model accounts for the optical properties of the tissue and also the spectroscopic properties of fluorescent dyes. The model is qualitatively validated against the the profiles of the spectra of multiple synthetic fluorescent dyes.
We are currently extending this pipeline to visualize the simulation data of in silico VSDI experiments to validate the simulation of the cortical activity for a large meso-scale circuit and also to visualize neuroglial cells. We are also working on accelerating the performance of the rendering workflow by providing a high performance distributed solution on multi-GPU visualization clusters based on the framework presented by Eilemann et al. [58].