Nuclear magnetic resonance (NMR) has proven of enormous value in the investigation of porous media. Its use allows studying pore-size distributions, tortuosity, and permeability as a function of the NMR sequence parameters [1]. This information plays an important role for characterising white matter (WM) tissue in vivo and non-invasively. A complete NMR analysis in silico involves the solution of the Bloch-Torrey equation over realistic domains. However, analytically solving this equation becomes intractable for all but the simplest geometries [1, 2]. To solve this limitation, numerical algorithms (e.g. [3]) are used for obtaining simulated signals in arbitrary domains. Nevertheless, geometrical models of WM are usually too simplistic for representing the real nature of tissue microstructure (e.g. parallel cylinders), raising doubts on the applicability of the results in live tissue.Develop statistical computational models of tissue blocks of WM microstructure based on histomorphometric data.

We based our analysis in axonal structures only, which were considered as generalised cylinders, i.e. non-straight cylinders with arbitrary cross-sections [4]. To do so, we create random instances of their axes by means of a random walk algorithm based on a multivariate von Mises-Fisher distribution [5]. This allows to consider a global axonal direction (i.e. the direction from one end of the axon to the other) and a local axonal dispersion (needed for representing local changes in the direction, i.e. tortuosity). We model the global axonal direction (i.e. the overall orientation considering all the axons passing through the voxel) using either parametric (e.g. Watson) or an arbitrary orientation density function (ODF), the latter described by the corresponding spherical harmonic decomposition [8]. Such representation is flexible enough to account for voxels belonging to different parts of the WM. Once the axes are defined, we assign a mean radii value for each axon according to a gamma distribution [6].

Then, we model each axon as a B-spline generalised cylinder using swept surfaces over the simulated axes [4]. This is a general technique providing flexibility for further generation of surface and volumetric meshes needed for numerical simulations (e.g. [3, 7]). Once the cylinders are built, we pack them avoiding intersections by random placement. The procedure is repeated until reaching a specific volume ratio. All computations were performed in MATLAB 2015a programming environment (The Mathworks, Natick, MA).

We tested the algorithm in different scenarios represented by the aforementioned probability distributions, requiring 3 minutes in average to obtain final results. Figure 1 shows two examples considering the same parameters but different global orientation density functions. In Figure 1a we show a random trial with global orientation represented by an arbitrary ODF (i.e. sampled from an arbitrary voxel extracted from the IIT atlas [8]), whereas in Figure 1b we considered a uniform distribution on the unitary sphere. It is noticeable the flexibility of the proposed technique for representing completely different axon dispositions.We developed a method for generating random samples of WM tissue from statistical data. The main advantage over existing methods is its ability to generate random samples of WM voxels based on histomorphometric data extracted from literature. This is a major step forward in representing axonal structures realistically, which is crucial for testing state-of-the-art diffusion MRI inverse models without being simplistic. Further work will be focused on increasing the complexity of the models to consider other neurites and intra-axonal structures. These models will then be used for simulating diffusion MRI signals and study its sensitivity to WM changes.