An Analytical Model for Estimating Water Exchange Rate in White Matter Using Diffusion MRI

Substantial effort is being expended on using micro-structural modeling of the white matter, with the goal of relating diffusion weighted magnetic resonance imaging (DWMRI) to the underlying structure of the tissue, such as axonal density. However, one of the important parameters affecting diffusion is the water exchange rate between the intra- and extra-axonal space, which has not been fully investigated and is a crucial marker of brain injury such as multiple sclerosis (MS), stroke, and traumatic brain injury (TBI). To our knowledge, there is no diffusion analytical model which includes the Water eXchange Rate (WXR) without the requirement of short gradient pulse (SGP) approximation. We therefore propose a new analytical model by deriving the diffusion signal for a permeable cylinder, assuming a clinically feasible pulse gradient spin echo (PGSE) sequence. Simulations based on Markov Random Walk confirm that the exchange parameter included in our model has a linear correlation (R2>0.88) with the actual WXR. Moreover, increasing WXR causes the estimated values of diameter and volume fraction of the cylinders to increase and decrease, respectively, which is consistent with our findings from histology measurements in tissues near TBI regions. This model was also applied to the diffusion signal acquired from ex vivo brains of 14 male (10 TBI and 4 normal) rats using hybrid diffusion imaging. The estimated values of axon diameter and axonal volume fraction are in agreement with their corresponding histological measurements in normal brains, with 0.96 intra-class correlation coefficient value resulting from consistency analysis. Moreover, a significant increase (p = 0.001) in WXR and diameter and decrease in axonal volume fraction in the TBI boundary were detected in the TBI rats compared with the normal rats.


Introduction
Diffusion weighted imaging (DWI) is an interesting non-invasive tool for measuring micro-structural and physiological parameters of tissue. However, because of several limitations and restrictions in both imaging and methodology techniques, the task of relating specific microstructures to the diffusion signal is not straightforward and requires many simplifying assumptions. On the other hand, for measuring and estimating micro-structural and physiological parameters of the tissue, such as cell size and water exchange rate between intra-and extra-cellular spaces, high bvalue diffusion imaging is generally required [1].
In contrast, there are also some high b-value diffusion approaches in which a model is constructed based on the tissue structure and physical rules such as Fick's laws to simulate diffusion and investigate its behavior under different conditions. Solving Fick's laws even in a simple tissue model requires extensive calculations. In these cases, Fick's law is solved using numerical techniques such as finite element [6], finite difference [7,8], or partial differential equations (PDE) [9]. For complicated structures, probabilistic methods like Markov Random Walks (MRW) may be used [10]. After simulating the water diffusion in the structure, it is possible to calculate the diffusion signal using proper mathematical relations. Thus, these models are good validation tools in building a physical model by ensuring that each of its parameters has an obvious effect on the diffusion signal [11,12]. Moreover, they can be used to estimate unknown parameters of the tissue. For example, the conventional DTI model combined with a physical simulation model of axons [13] is used to estimate white matter tissue parameters such as density and myelin sheath thickness [14].
Other approaches utilize and optimize the parameters of the MRI pulse sequence in order to make it sensitive to specific microstructural parameters [15]. For example, double pulsed-field gradient (d-PFG) MRI scheme is used to measure pore diameter [16] and double-wave-vector imaging protocol is used to estimate the density and size of cells even in irregular distributions [17]. The phase of the diffusion signal also contains useful information about the tissue such as fiber orientation [18]. On the other hand, high field diffusion imaging is also an interesting tool to measure cell size because of the diffraction patterns seen in the diffusion signal in high field gradients [19]. This approach is employed to measure the size of microstructures in a porous medium [20,21]. However, it hardly can be applied to white matter fibers due to its limitations in both imaging and modeling methodologies [21]. Another example of high field imaging application is the work of Jespersen et al. [22] in which a two-compartment diffusion model is proposed to measure the density of neurites (axons and dendrites) for a normal rat.
Presently, most researchers use muli-compartment models for modeling of diffusion signal acqired by MRI [1,[23][24][25]. Typically, they consist of at least two compartments: one compartment is devoted to diffusion in the extra-cellular or extra-axonal space, known as hindered diffusion or fast diffusion. The other part is diffusion in the intra-cellular or axonal space, known as restricted diffusion or slow diffusion. Different approaches define and use different parameters and different simplifying assumptions. The CHARMED (Composite Hindered And Restricted ModEl of Diffusion) introduced by Assaf et al. [1] is a well known model that has been used in several studies to measure tissue parameters, especially axon diameter. For example, Assaf et al. [23] estimated axon diameter distribution of freshly excised porcine optic and sciatic nerves using this model and high b-value MRI perpendicular to the axons. However, one of the first attempts to measure axon diameter distribution in vivo was conducted by Barazany et al. [25] for the corpus callosum fiber tracts of rats by combining DW imaging and a new analytical method called AxCaliber proposed by Assaf et al. [24]. This work is an extension of the CHARMED model which has an additional compartment for isotropic diffusion, accounting for cerebrospinal fluid (CSF). In this approach, the diffusion gradients also need to be perpendicular to the fiber orientation. To reduce the dependence on the fiber orientation, Alexander et al. [26] proposed a four compartment model consisting of intra-axonal space, extra-axonal space, CSF, and glial cells combined with a multi-shell HARDI data acquisition. The scalars defined in this method for estimating the axon diameter and density are invariant to the orientation of the axons. Another limitation of this kind of modeling is that a single orientation for the axons is assumed which limits the analysis to the voxels belonging to single fibers. Therefore, Zhang et al. [27], using a simplified version of the CHARMED model, brought in the effect of orientation dispersion of the axons to the modeling by assuming Watson distribution for the axon's orientation. Furthermore, by modifying their model they [28] proposed a three comportment model combined with a clinically feasible two-shell HARDI protocol to estimate the neurite orientation dispersion and density (NODDI). Hence, this approach can be applied to a wide range of voxels within the brain. As a final point here, Panagiotaki et al. [29] did an experimental study to investigate the performance of 47 multi-compartment models of the white matter with no exchange. They concluded that by combining a diffusion tensor (DT), a cylindrical, and a spherical component, the constructed model describes the data more accurately. Moreover, considering a distribution rather than a single value for the axon's radius decreases the stability of the fitting process.
It should be mentioned here that Ball and stick [30] is another kind of modeling which consists of two terms: one for modeling the isotropic diffusion, ball, and the other is accounting for the fibers, zero-radius stick(s). Although this kind of modeling looks similar to two compartment models like CHARMED, it is mainly useful for resolving the fiber orientation in crossing regions.
An interesting parameter of the white matter is the water exchange rate between the intra-and extra-axonal spaces, which also has an effect on the diffusion signal [31]. Karger et al. [32] added the exchange effect into the NMR magnetization differential equations to derive proper relations to measure exchange times. One of the first analytical models including exchange rate was introduced by [33] for the bovine optic nerve using Karger relations. There are also other studies which investigated the effect of the exchange rate based on the Karger framework [12,31,34,35]. Injuries or abnormalities of the cerebral tissues can change the exchange rate [35]. Moreover, it is shown using Monte Carlo simulation that in a two-compartment diffusion model, the exchange rate through a permeable membrane should be included so that the estimated parameters are close to their actual values [36]. Usually, these models are valid when the diffusion gradient is applied perpendicular to the tract. The short gradient pulse (SGP) condition also needs to be valid when applying the NMR imaging protocol. On the other hand, some techniques modify and optimize the acquisition pulse sequence to make it sensitive to the exchange rate, such as double-PGSE diffusion exchange spectroscopy [37], filter exchange spectroscopy (FEXSY) [38], and apparent exchange rate (AXR) imaging [39].
Another way to import the exchange effect into the modeling of the diffusion signal is to solve Fick's Laws for a typical structure with proper boundary conditions, which considers permeable membranes [40]. However, even for an impermeable boundary, solving Fick's equations is difficult and needs simplifying assumptions such as SGP or identical diffusion time, D, and diffusion pulse width, d, [41] which limit the applicability of the methodology for a typical human MRI scanner.
Therefore, the aim of this work is to present a new analytical model for the diffusion signal which includes axon diameter and axonal volume fraction as well as water exchange rate. The diffusion signal is calculated mathematically for a pack of parallel cylinders with permeable membranes, assuming a practical PGSE sequence. The model has been validated and confirmed using artificial data with Markov Random Walk simulation, DWMRI, and histological data from the ex vivo rat brain. The model is also used to compare TBI rats with normal ones with respect to the estimated parameters, especially in the corpus callosum.

Materials and Methods
In this section, we derive the diffusion signal assuming a PGSE sequence perpendicular to a typical cylinder having permeable membrane. Then the other parts of the model, accounting for the diffusion signal of parallel and extra axonal spaces, are combined accordingly. After that, the methodology for evaluating the proposed model using artificial data simulated by Markov Random Walk is described. Finally, the DWI protocol and the staining process employed in our experiments are explained.

a) Background and Motivation
The general equation for relating the diffusion signal to the movement of the spins during diffusion time, D, is in which (w 2 2w 1 ) is the phase change due to the spin displacement along the direction of the diffusion gradient. To solve equation (1), the probability function of the phase change has to be known. The simplest way is to assume a Gaussian distribution for this function; therefore we will have This assumption is valid when diffusion is un-bounded or free in the underlying medium like inside the extra-cellular space. However, usually the researchers use the above equation also for bounded media like the axons in order to simplify the problem [41]. Under this assumption, the problem is reduced to calculate the ''mean square phase change due to the spin displacement'', i.e., in which c is the gyromagnetic ratio, d and D are the pulse width and diffusion time, respectively. Moreover, G is the amplitude of the diffusion gradient, which we assumed along the x-axis here. This also requires the probability function (propagator) of spin displacement, P(r|r'), which comes from solving Fick's second law under proper initial and boundary conditions consistent with the underlying structure.
Another way to solve equation (1) is to assume d,,D (short gradient pulse, SGP, approximation) and then use the Fourier representation of the problem [42]. Therefore, the Gaussian assumption for phase change does not have to be imposed. However, the propagator function P(r|r') is still required.
SGP approximation can hardly be achieved in existing clinical and even research MRI scanners. On the other hand, the first method imposes the Gaussian distribution on the model, which restricts the number of unknown parameters. For example, Neuman [41] calculated the perpendicular part of intra-axonal diffusion signal as where q = cGd. Here, just one unknown parameter can be estimated if the images are acquired using a fixed value of D. For instance, in the CHARMED model [1], only the radius, R, is estimated using this part, while D in is estimated using the parallel part of the signal. To add another parameter to the model, such as exchange rate, there are two options; either we can define and use a new imaging protocol with different diffusion times, like [37], or we can build our model in a way that does not require Gaussian assumption for phase change and also does not require SGP approximation. Although the second approach is more realistic, to our knowledge, there is no calculation trying to use it because it necessitates non-solvable integrals. Here, we built an analytical model including a parameter accounting for the exchange rate using proper mathematical relations without assuming Gaussian distribution for the phase change or the SGP approximation.

b) Proposed Model
The phase change can be approximated by the summation instead of the integral: where N..1. Assuming diffusion gradient is along the x-axis, the more general form will be Dw~2 pq q N : where,ṽ v~2 pq q N ,r r D (k)~r r(Dz(k{0:5) d N ), and r r(k)~r r((k{0:5) d N ). Then, using equation (1), the diffusion signal is calculated by the following integral: ::: ð r 1 dr 1 dQ 1 :::r N dr N dQ N :r D1 dr D1 dQ D1 ::: |P(r r 1 ,:::,r r N ,r r D1 ,::: in which dr r~rdrdQ. Since the diffusion is a Markov random process, each movement of a spin depends only on its starting position. Therefore, after replacing dr 1 dQ 1 :::dr N dQ N by dRdW and expandingP(r r 1 ,:::,r r N ,r r D1 ,:::r r DN ), equation (7) is simplified to: ::: Pr r 1 |P(r r 1 Dr r 2 ):::P(r r N{1 Dr r N )|P(r r N Dr r D1 )::: Note that the imaginary part of equation (7) is equal to zero. Now we need the probability function (propagator) of spin displacement, P(r 0 |r). This function has been calculated for planar, spherical, and cylindrical boundary conditions [41,42] and used in several calculations. However, few of them consider the effect of permeable membrane in the boundary condition [40]. For a cylinder of radius R, the 1-D propagator function perpendicular to the axon is [42]: which is the solution of the Fick's Law, subject to the following initial and partially permeable membrane boundary conditions: where D in is the diffusion coefficient of the intra-axonal space and M(m/s) is the permeability coefficient. J n (.) is the n th order Bessel function and a nm is the m th positive root of the following equation which comes from combining equations (10) and (12).
The right hand side of the above equation is a dimensionless quantity which we call absorbing coefficient, h = R.M/D in . This quantity is directly related to the exchange rate (1/s), defined as the fraction of spins exchanged between intra-and extra-axonal space through the membrane per second, because each spin absorbed by the membrane moves to the extra-axonal space and there will be also a corresponding spin which is transported into the intra-axonal space in order to preserve the equilibrium. The coefficients B nm can be derived using the orthogonality property of the Bessel Functions and the above boundary condition: By solving equations (8), (9), (13) and (14), the intra-axonal compartment of the diffusion signal perpendicular to the axon is achieved. Starting from equation (8), by taking the integrals over angels, Q, using the following equation the equation (8) is reduced to: E(q,D)~ð ::: Pr r 1 |P(r r 1 jr r 2 ):::P(r r N{1 jr r N )|P(r r N jr r D1 )::: in which Pr r 1~1 pR 2 and the other probability functions are derived from equation (9). The most dominant term in that equation is the one having a 01 [40,43], therefore, P(r r k Dr r kz1 ) d=N , for k~1:::N{1 where and the simplified boundary condition will be After inserting (17) and (18) into (16), we have Unfortunately, only the second integral has a closed form which is To get a closed form for the first integral, we assume thatJ 0 (vr)%1{ v 2 r 2 4 which is valid because N is considered large enough in our analysis. Using this approximation and also the following integral equation, the closed form of the first integral appeared in equation (21) will be R 2 Finally, by simplifying equation (21) using equations (20), (22), and (24) we have Therefore, the final form of the signal which we use in our model is: where a 01 and h are related by equation (23).
Here we should note again that in equation (9) we omit the terms with n.0 and m.1 from the summation as also assumed by [40,43]. From equation (13) it can be shown that a nm R$3.83 for all values of n and m except (n,m) = (0,1). Therefore, to ensure that this approximation is valid, we should have: Thus, there is an upper limit for N, which is Nv3:2 D in d R 2 . Combining this limit with the lower limit introduced in equation (5), we have: For example, since in our experiments, the typical values of D in , R, and d are 0.25 mm 2 /ms, 0.5 mm, and 10 ms, respectively, a reasonable value of N is 30 which is chosen in our analysis.
For special cases when the SGP approximation is valid, i.e., d,,1, we can set N equal to 1. Then, equation (26) is reduced to the following form of the diffusion signal, which is similar to the form used by [40]: Moreover, if the absorbing effect, h, (or the exchange phenomenon) is neglected from the equation, we have: which is consistent with the form of the signal introduced by [43], for long D. Similar to [1], if the diffusion gradient, G, is decomposed into parallel and transverse directions of the axon, i.e., G || and G H , and since q = cdG, the attenuated signals due to diffusion in and out of the axons are as follows assuming Gaussian profile for E in || , E out

||
, and E outH . The overall diffusion signal is a weighted summation of both compartments based on their volume ratios: in which v is the intra-axonal volume fraction ratio. Also, if the Johnson noise model is considered for the measurements [1], the final model for the measured diffusion signal is.
In summary, our two compartment model consists of four parts to model the diffusion signal. Among these parts, the effect of membrane permeability is modeled explicitly only by E inH which takes into account the spins inside of the axon absorbed by the membrane. This means that those exchanging spins which leave or enter the intra-axonal space are modeled by the other parts which assumed to have Gaussian distribution. As a result, this assumption may affect the estimation of the other parameters of the model especially the apparent diffusion coefficients as in current nonexchanging two compartment models. However, as stated before, the absorbing coefficient, h, included in the proposed model is a quantity which can be used to estimate the exchange rate.

c) Model Fitting
Having measurements and the proposed model, the remaining task is to fit the model to the data. In this work, we use a non-linear least square fitting procedure to estimate unknown parameters (D in , D out,|| , D out,H , R, h, v, and g). For the fiber direction, we use the principal eigen-vector estimated from the DTI model. Note that since the proposed model is built for a specific configuration of fibers (parallel axons with nearly identical diameters), it is valid just for a special group of voxels. Therefore, we apply the model to the voxels having high FA values which are more likely to belong to a single fiber tract. We also assume that in a single fiber, the axons are parallel and nearly the same size. The fitting procedure is implemented in MATLAB (MathWorks, Natick, Massachusetts, USA) by employing its optimization toolbox. A block diagram of the proposed model for driving the diffusion signal is shown in Figure 1.

d) Markov Random Walk Simulation
Validation of a model can be conducted using simulation data. In this work, we use Markov Random Walk (MRW) to simulate the diffusion process and create an artificial diffusion signal in a voxel. An advantage of MRW is that it can be applied even to complex configurations of the fibers. Nevertheless, a disadvantage of MRW is its huge computational complexity. The aim of this simulation is to solve the diffusion equation for a fiber.
Consider a pack of parallel cylinders of equal radius, R, (rolling as axons) distributed in a regular manner (see Figure 2(a)). Thus, the medium is partitioned into two components: intra-and extraaxonal spaces having the diffusion coefficients of D in and D out , respectively. To simulate diffusion, first we construct a 3D grid of cells to discretize the space. In this work, we assume a voxel with dimension of 24624648 mm discretized (separated) into 806806160 cells (Dx = Dy = Dz = 0.3 mm). Axons are considered as parallel cylinders in the z-direction with a radius of 1.5 mm. Axonal volume fraction, v, is defined as the ratio of the intraaxonal space to the total volume. In this work, we assume two configurations with v = 0.27 and v = 0.46.
To simulate the random walk, each spin starts from its initial position and moves with a step size, dr = !(D.dt/6), in one of the 6 directions along the Cartesian axes, i.e., 2x, x, 2y, y, 2z, z, randomly. Note that the step size, dr, is different from the cell size, Dx, Dy, or Dz. So, based on its position and the diffusion gradient's direction and magnitude, the phase of that spin is updated by.

w(tzdt)~w(t)zcg g(t):r r(tzdt) ð34Þ
with the initial value of Q(t = 0) = 0 and time step, dt, equal to 2 ms. Then, the attenuated signal due to diffusion at echo time, t = TE, is calculated using the following relation: At t = 0, 10 6 spins are distributed uniformly on a disk which is placed in the center of the voxel (see Figure 2(b)). The radius of the disk is selected such that the ratio of the initial spins inside the cylinder to the outside ones is equal to the intra-to extra-axonal volume ratio. We should emphasis that it is not necessary to initialize the spins across the whole substrate because in our simulation the configuration is symmetric. Moreover, in this way, the likelihood that a spin reaches to the boundary of the voxel is significantly decreased which is desirable in the numerical calculations. To calculate the diffusion signal, we consider a typical gradient spin echo sequence as depicted in Figure 3. The sequence parameters are chosen similar to our ex-vivo animal experiment with d = 10 ms, D = 18 ms, and q = 0.025, 0.05, 0.075, 0.1, 0.125 1/mm, corresponding to 357, 1437, 3240, 5767, 9016 s/ mm 2 b-values, respectively. However, the number of diffusion gradients is set to 6 for all five shells to reduce computation time.
To include the exchange process into the simulation, the spins reaching the boundary may pass it, based on a jump probability, s, which is calculated from the following equation, in which P m is the membrane permeability, Dm is the membrane thickness, and k represents the relative concentration of the inside/ outside of the axon to the free water which we set it to 1 in this Water Exchange Rate (WXR), 1/sec, is calculated as the ratio of the spins passed the boundary to the total number of initial spins per second. Moreover, to study the performance of the proposed model more realistically, white Gaussian noise is added to the signals (SNR = 40) and each investigation is repeated 50 times. In Figure 4, the final distribution of the particles is shown in a cross section of the structure. In this simulation, the permeability is equal to zero; therefore, the particles did not pass the boundaries. In Figure 5, the mean and STD values for the estimated parameters h/D, r9, v9, and the diffusion coefficients are shown versus WXR. Note that h/D is the absorbing coefficient normalized to the diffusion time. It can be seen that in Figure 5 (a) h/D has a significant (p = 0.0001) linear correlation to the WXR, with R 2 -values of 0.91 and 0.88 for v = 0.46 and v = 0.27, respectively. Also, the intra-class correlation analysis within 95% confidence interval shows 0.968 and 0.967 agreement (1:1 match) (p = 0.0001) between h/D and WXR for these two configurations, respectively. On the other hand, the diffusion coefficients have also been well estimated using the model. However, for radius, r, and volume fraction, v, the estimated values change as the WXR increases. Note that v9 actually corresponds to the ratio of spins trapped inside axons, which decreases by increasing WXR. On the other hand, by decreasing v9, one may expect that the radius should decrease too. However, it should be noted that in our model, the number and geometrical configuration of the axons is not considered. Therefore, increasing the radius does not necessitate increasing the intra-axonal space. Moreover, one reason the estimated radius, r9, is increased by increasing WXR is that when more spins pass the axon's boundary, the effective apparent radius increases. It should be stated that compensating for these two effects is almost impossible using two-compartment models. However, it may be solved if the number and geometrical configuration of the axons are known and imposed into the model using proper relations.

e) Experimental Design
To estimate unknown parameters using the proposed model, there should be appropriate diffusion MRI data. First, data should be acquired using more than one b-value (multi-shell diffusion   imaging) so that the profile of diffusion in the radial directions can be revealed. Secondly, the diffusion time, D, should be large enough so that the spins have enough time to reach the boundaries.

All Experimental Procedures Were Approved by the Institutional Animal Care and Use Committee (IACUC) of Henry Ford Health System
Fourteen male Wistar rats (weight: 270-300 g) were divided into two groups, traumatic brain injury (TBI) (n = 10) and a normal group without neurological injuries (normal group, n = 4). All TBI rats were subjected to controlled cortical impact (CCI) [44] and sacrificed six weeks after TBI. The brain was perfused with heparinized saline [44]. Ex vivo MRI scans were performed seven weeks after TBI. Normal rats were sacrificed and followed the same procedures and measurements as TBI rats.
3) Histological staining. The procedures to prepare brain sections were the same as previously reported [44]. Bielshowski and Luxol fast blue staining were used to identify reticular fibers (i.e. neurofibrils and neurofibrillary tangles) and myelin, respectively. For immunohistochemical staining, slides were placed in 20% silver nitrate in the dark, and ammonium hydroxide was added until the tissues turned brown with a gold background, and then sodium thiosulfate was added. The slides were stained for Luxol fast blue, washed in 95% alcohol and placed in lithium carbonate. Nuclei are colorless; myelin is blue and axons are black.
4) Histological analysis. Two sections of each brain were selected for taking light microscopy images with different magnification lenses (2x and 100x). The images were taken from seven parts of the corpus callosum (CC) as shown in Figure 6. If the tissue related to any of these ROIs is corrupted due to TBI or cutting process, the nearest region belonging to CC is considered h/D has a significant (p = 000.1) linear correlation to WXR, with R 2 -values of 0.91 and 0.88 for v = 0.46 and v = 0.27, respectively. Also, the intra-class correlation analysis within the 95% confidence interval shows 0.968 and 0.967 agreement (1:1 match) (p = 0001) between h/D and WXR for these two configurations, respectively. It is also seen that the diameter increases while the axonal fraction decreases as the WXR increases. doi:10.1371/journal.pone.0095921.g005 for the photography. in Figure 7, typical images of ROI #5 for one of the TBI samples are shown.
The 100x magnified images were used for measurement of axon diameter and intra-axonal volume fraction. To measure the axon diameter, in each image at least 10 axons were selected visually and the number of pixels across the axons were counted. Since each pixel represents a 0.196 mm by 0.196 mm region in the actual tissue, axons thinner than this value may not be included in our calculations. To estimate the volume fraction, the images were converted to gray scale images. Then, after color inversion (to make the axons bright instead of dark) and histogram enhancement, the axons were extracted by thresholding the images using Otsu's method [46]. The soma and nuclei were also eliminated from the final binary image using morphological operations.

f) Validation of the Proposed Model Using Histological and DWMRI Data
As we described in Section (b), the brains were scanned before staining using the aforementioned protocol. The resulting images were smoothed using a 363 Gaussian filter to reduce the noise. For each brain, the corpus callosum (CC) fiber tract was segmented manually and divided into 3 parts: left, center, and right for comparison (see Figure 8). The remaining task for model validation is to find the corresponding regions between the MRI and histological images. This can be done by registration of each histological section to its corresponding T2-weighted image.
However, even using a nonlinear registration algorithm, these two images will not be matched perfectly due to the significant intrinsic differences between these two modalities. For example, the MR imaging slice thickness is 1 mm while the histological section thickness is 6 mm. Moreover, there are some distortions and ruptures in the sections caused by the staining process. Therefore, a reliable way to find the corresponding regions is to do it visually. In our validation step, the center part of the corpus callosum of the normal cases is considered because this part is thick enough to assume that it is almost the same in the range of 1 mm. Also, correspondence of this region between MRI and histology can be established more reliably than the other regions. Needless to say, we can only validate the model with respect to the estimated structural parameters, i.e., the diameter and axonal volume fraction.

A) Histological Measures
The measured values for diameter and axonal volume fraction of each ROI in the histological section for each case are shown in Tables 1 and 2, respectively. The corresponding plots of the diameter and axonal volume fraction versus ROI with error-bar are also shown in Figure 9. As seen, there is a significant difference between the normal and TBI rats in ROI#5 close to the TBI    boundary, both for diameter (p = 5.7610 213 ) and axonal volume fraction (p = 0.0033). In other words, the injury in ROI#5 after TBI has increased axonal diameters and reduced axonal volume fraction compared to the uninjured animals. As discussed in Section (d) of 'Material and Methods', these findings are consistent with the simulation results of our proposed model, where we found that by increasing WXR, the axonal volume fraction decreases while the axon diameter increases. Therefore, by considering the simulation results and histology findings, it can be claimed that the WXR increases in the injured brain. Table 3 shows the estimated values of the diameter, 2R, and axonal volume fraction, v, for the center part of the CC of normal cases using our proposed model. By comparing this table to the  histological values shown in Tables 1 and 2, especially ROI#4, it can be seen that v is underestimated (about 30%) while R is overestimated (about 20%). However, the results are consistent because the consistency analysis between the results of histology and the proposed model reveals that the intra-class correlation coefficients within 95% of confidence interval are 0.970 (p = 0.008) The observed differences between the results of the model and the histological measurements can be caused by several factors. These factors generally originated from the limitations in modeling (simplifying assumptions like parallel, same size axons), imaging (like noise), and also histology staining process. The latter factor includes the inevitable unwanted incidences which may occur during cutting and staining such as corruption, contraction, or expansion of the tissue. On the other hand, measurement of the radius and axonal volume fraction has its own limitations. For instance, the pixel size on the 100x-magnified images is 0.196 mm which imposes an uncertainty of about 0.2 mm to the values of measured axon diameter. Moreover, the thickness of the sliced sections (which is 6 mm) has not been accounted for measuring the volume fraction.

C) Application of the Proposed Model to Discriminate TBI Tissues
To show the capability of the model to discriminate between the injured and normal tissues, three parts of the CC were segmented manually as depicted in Figure 8 and investigated using the proposed model. In Figure 10, the error-bar diagrams (mean with STD values) of the estimated parameters, h/D, 2R, and v in different ROIs are plotted for all cases including normal and TBI rats. There is a clear contrast with p-value of 0.001 between the injured (right) and normal (left) parts of the CC, especially for the map of h/D parameter (related to the water exchange rate, WXR) which significantly increases in the TBI boundary compared with that in the controls. Also, as expected, the diameter increases and the axonal volume fraction decreases in the injured fibers, which is consistent with our simulation results.
The proposed model can also be employed to investigate other diseases affecting the white matter fiber tracts such as multiple sclerosis (MS). Especially, the estimated values of h in the fiber tracts close to or within the damaged area may be used as a quantitative measure to show the progress of the disease. For example, this parameter may demonstrate the level of demyelination since there is a meaningful relation between the myelin sheath and the exchange rate. Therefore, as the proposed model can be applied to the data acquired using clinical scanners, our model provides a new insight that may be relevant to clinical diagnosis and treatment of neurological diseases.
Although the methods presented in this work are based on high b-value DWMRI, it is worthy to discuss the effect of changing micro-structural parameters on the DTI parameters such as fractional anisotropy (FA) and mean diffusivity (MD). Generally, increasing WXR may decrease FA while increase MD. On the other hand, if the intra-axonal space increases, the diffusion inside the axons becomes less 'restricted', however, the diffusion outside of the axons becomes more 'hindered'. Therefore, if the hindered part is the dominant part of the diffusion signal (which is true in common DTI protocols), by increasing the intra-axonal space the FA is increased while the MD is decreased, and vice versa. It should be implied that the true effect of variation of the axonal space on the FA and MD values is dependent on several factors such as imaging parameters and tissue specifications.
Brain injury, such as TBI, remains a leading cause of mortality and disability among children and young adults. Current research in brain injury has been restricted to acute neuroprotection  treatment with a short treatment window [47][48][49]. In these investigations, functional recovery after TBI may be driven by neuronal remodeling. Currently, investigation of neuronal remodeling after brain injury has been dominated by traditional DTI such as FA and fiber tracking [50][51][52][53][54][55]. However, conventional DTI produces an anomalous result, showing an overall lowering of FA despite the presence of highly-oriented tissue [3,44,[56][57][58][59][60].
Currently, there are few published studies to target the issue related to low FA in areas with highly-oriented tissue [44,57,58,60,61] and no investigation has been published in white matter remodeling after brain injury, such as TBI, using MRI axonal permeability and diameter. Our data demonstrated that physiological measurable parameters, MRI axonal permeability, diameter, and density are sensitive in detecting axonal changes after TBI and could be very useful in the evaluation of axonal damage and remodeling after neurological diseases, such as TBI, stroke, and MS.

Conclusion
In this work, a new analytical diffusion model was proposed for estimating micro-structural (axon diameter and volume fraction) and diffusivity (diffusion coefficients as well as water exchange rate) parameters of white matter fiber tracts using DWMRI. Therefore, one of the specifications of the proposed model is that it includes the absorbing coefficient of the intra-axonal space. This coefficient is linearly related to the water exchange rate (WXR), which has not been fully investigated in previous calculations. Moreover, the proposed model does not require short gradient pulse (SGP) approximation and thus can be applied to the data acquired using common PGSE sequences in a clinical scanner. The simulation results using Markov Random Walk showed that the included parameter, h, in the model has a significant, proportional relation to the actual water exchange rate, WXR. Moreover, the simulations revealed that by increasing WXR the axon volume fraction decreases while the diameter increases. Consistent with our model, histological analysis of TBI brain rats revealed that in the injured tissues, the volume fraction reduced and the axon diameter increased. Therefore, according to these findings, it can be concluded that the WXR increases in injured tissues, which is reasonable, because the cell membrane, especially the myelin sheath, is damaged. The estimated values by applying our proposed model to the MRI data of normal rats were relatively comparable to their corresponding histological measurements. Nevertheless, the observed differences can be caused by several factors originating from the limitations in modeling, imaging, and also the histology experiment. As shown in this work, the proposed model can be used to discriminate TBI and normal tissues with respect to the estimated parameters, especially WXR. One of the limitations of the proposed model, as usually exists in other similar models, is that it can only be applied to specific regions of the white matter where the axons are parallel and the same size. This issue can be addressed in the future by adding more compartments to the model as well as improving the MRI data acquisition.