With the proposed model (Figure 1), both the MRI and EEG data are described as a function of the distribution of axonal radius [P(r)] and of the axonal g-ratio [g(r)] within a given white matter tract. The axonal radius distribution is assumed to be a Gamma distribution (Sepehrband et al., 2016):

where r is the axonal radius. M represents the mode, i.e., the peak of the axonal radius distribution and the parameter θ represents the width of the right tail of the axonal radius distribution, a measure of the number of large axons in a white matter fiber tract.

Histological studies have shown that larger axons exhibit comparatively thinner myelin sheaths, i.e., larger g-ratios (Ikeda and Oka, 2012; Gibson et al., 2014). From data obtained in the peripheral nervous system of the rat, the radius dependence of the g-ratio was shown to follow (Ikeda and Oka, 2012): g_REF(r) = 0.22 log(2r) + 0.508. However, to facilitate the mathematical manipulation of the biophysical model (Eqs. 3 and 5 below), we write the radius dependence of the axonal g-ratio as:

With this power law, the exponent α represents the slope of the radius dependence of the axonal g-ratio, while the parameter β is a scaling factor. Unlike the mathematical expression for g_REF, this power law does not include an offset term. To verify its validity, we fitted this power law with the reference relationship (g_REF) of Ikeda and Oka (2012). The result shows an excellent level of agreement between both expressions (R² = 0.99), with α = 0.18 and β = 0.57 (Figure 2A). We also verified the applicability of the power law in the human central nervous system, where the higher axonal g-ratio requires the addition of a systematic offset to g_REF. We estimated this offset to be 0.14, assuming a g-ratio of 0.7 for an axonal radius of 0.9 μm, based on studies showing g-ratio estimates between 0.65 and 0.79 (Mohammadi et al., 2015; Stikov et al., 2015) across the range of axonal radius observed in the human brain (Aboitiz et al., 1992; Caminiti et al., 2009; Liewald et al., 2014). Fitting of the power law (Eq. 2) with the reference relationship (g_REF) after addition of this offset leads to an excellent agreement (R² = 0.99), with α = 0.14 and β = 0.71 (Figure 2B). The latter values will be used when the parameters α and β are set constant to allow for the estimation of other model parameters (see section “Model Parameters”).

In agreement with West et al. (2016), the MRI-measured g-ratio is written as an ensemble average of the axonal g-ratios within each image voxel, weighted by the axons’ cross-sectional area:

where R is the fiber radius [R = r/g(r)]. From Eqs. (1) and (2), Eq. (3) becomes (see Supplementary Appendix A):

As described in Waxman and Bennett (1972), axonal conduction velocity (v) can be derived from the morphological properties of axons using: v[m/s]=pd[μm]g, where p (∼5.5–6.0) represents the contribution of additional axonal factors to the propagation of action potentials (e.g., length of Ranvier nodes, electrical properties of the myelin membranes). Assuming an equal contribution from all axons to the conduction velocity V measured with EEG, we obtain:

From Eqs. (1) and (2), Eq. (5) becomes (see Supplementary Appendix A):

The proposed model (Eqs. 4 and 6) includes 4 free parameters, pertaining to the distribution of axonal radius (M and θ) and to the radius-dependent axonal g-ratio (α and β). However, our proposed model uses only two data types acquired in vivo (MRI-measured g-ratio and EEG-based axonal conduction velocity). It is therefore necessary to set two parameters to reference values.

In general terms, the choice of model parameters to estimate should take into consideration the neural mechanisms underlying each application study of this model. The exponent α of the power law g(r) represents the rate of change in thickness of the myelin sheath with axonal radius and may be a parameter of interest when considering changes in neuronal shape that differentially affect axons of different sizes. In contrast, the scale parameter β equally affects axons of all sizes. Concerning the distribution of axonal radius, we highlight that histological studies across white matter tracts and animal species have reported that, for reasons pertaining to brain size limitations and metabolism, the mode M of the axonal radius distribution remains largely constant (Tomasi et al., 2012; Liewald et al., 2014). This motivates setting M to a constant value from the histological literature and estimating the tail parameter θ from the in vivo data.

The choice of constant model parameters may also be guided by the impact of inaccurate constant values on the estimated morphological features. The bias of the parameter estimates arising from inaccuracies of 10% in α, M or β was evaluated using numerical simulations (with the procedure described in methods section “Estimation of Axonal Morphology From in vivo Data”). Across a range of plausible g_MRI and V, an inaccuracy of 10% in α leads to an average bias of ∼10 and ∼15% on the estimated M and θ, respectively (Supplementary Figure 1A). Similarly, an inaccuracy of 10% in M leads to an average bias of ∼1 and ∼22% for β and θ, respectively (Supplementary Figure 1B). However, an inaccuracy of 10% in β leads to an average bias of ∼155% and ∼92% for M and θ, respectively (Supplementary Figure 1C).

From the biological and numerical considerations above, in this study, we chose to set α to 0.14 (see section “Axon Morphological Properties”) and M to 0.40 μm (Tomasi et al., 2012; Liewald et al., 2014) and to estimate the parameters β and θ.

Numerical simulations were conducted using Eqs. (4) and (6), to examine the values of the model parameters (M, θ, α, and β) obtained from combinations of the in vivo data (g_MRI and V) using the procedure described in methods section “Estimation of Axon Morphology From in vivo Data”, and vice versa.

In order to highlight the range of in vivo data values compatible with the proposed model, estimates of M and θ were computed across a large range of g_MRI and V values, setting α = 0.14 and β = 0.71 (section “Axon Morphological Properties”). The range of in vivo data values compatible with the proposed model was determined by comparison of the M and θ estimates against reference literature values.

Subsequent numerical simulations were conducted by estimating the parameters β and θ across the range of compatible in vivo data values, setting α = 0.14 and M to 0.40 μm (section “Model Parameters”). In particular, we examined the range of the parameters β and θ across values of g_MRI and V reported in the literature.

The variability of the model parameter estimates in the presence of noise was also investigated. Simulated estimates of the in vivo data g_MRI and V were computed from combinations of θ and β values using Eqs. (4) and (6). Noise was added to the computed g_MRI and V values with a standard deviation of 0.03 for g_MRI and 0.50 m/s for V, representative of intra-subject variability in in vivo data. To replicate in vivo conditions, 700 samples of g_MRI and one sample of V were taken from the resulting distributions of g_MRI and V, and estimates of θ and β from noisy data were calculated. This process was repeated 2,000 times and the standard deviation of the θ and β estimates across repetitions was computed as a measure of their variability.

Finally, we set out to investigate the effect of using a group averaged conduction velocity rather than subject-specific velocities. We selected 15 samples of g_MRI and V values from distributions with means of 0.70 and 10 m/s and standard deviations of 0.05 and 0.80 m/s, respectively, representative of inter-subject variability in in vivo data. From these simulated in vivo data, reference values of β and θ were calculated. Estimates of β and θ were also computed from the average of the 15 samples of V. We then estimated the bias between the reference β and θ values and those obtained from the average value of V.

We acquired data from 17 right-handed healthy volunteers. All participants had normal or corrected-to-normal vision and hearing and had no history of psychiatric or neurological disorders. Participant handedness was evaluated with the Edinburgh Handedness Inventory (Oldfield, 1971). All participants gave written informed consent and received 80 Swiss Francs as monetary compensation. The study was approved by the local ethics committee.

MRI data quality was assessed using the Motion Degradation Index (MDI) described in Castella et al. (2018) and Lutti et al. (2022). MDI values were ≲ 4 s^–1 for the PD- and T1-weighted raw images, and ≲ 5 s^–1 MT-weighted raw images, indicative of good quality images (Lutti et al., 2022; see Supplementary Figure 2). Three participants were excluded due to artifacted EEG recordings. The final sample consisted of 14 participants (6 females; age = 27.14 ± 3.86 years). Three participants had left eye dominance, as determined by the eye viewing an object at a distance when the participant looks through a small opening (Miles, 1930). An overview of the data processing pipeline of the in vivo data is shown in Figure 1A.

Estimation of axonal morphological features from the in vivo MRI and EEG data was implemented using MATLAB-based custom-made analysis scripts. The MRI-measured g-ratio samples along the visual transcallosal tract and the estimate of the IHTT were used to estimate model parameter values using Eqs. (4) and (6) (see section “Axon Morphological Properties”, Figure 1B). This was achieved using MATLAB’s non-linear least-square routine (lsqnonl) with a trust-region-reflective minimization, which minimizes the sum of the squares of the residuals. The initial conditions used to ensure convergence of the fitting routine were set to β = 0.70 μm^−α and θ = 0.10 μm.

All codes are available on our online repository: https://github.com/LREN-physics/AxonalMorphology.

We investigated the range of in vivo data compatible with the proposed model by considering the values of the model parameters M and θ computed from each combination of g_MRI and V (Figure 3). Combinations of large values of g_MRI and low values of V lead to unrealistically low values of M (∼10^–14 μm). Conversely, small values of g_MRI and large values of V lead to excessively large values of M (∼1μm) and low values of θ (∼0.001 μm) (Aboitiz et al., 1992; Caminiti et al., 2009, 2013; Tomasi et al., 2012; Liewald et al., 2014).

Figure 4 shows the range of the model parameters θ and β across a large range of g_MRI and V values. According to the literature, frontal transcallosal white matter exhibits MRI-measured g-ratio values of ∼0.62 (Mohammadi et al., 2015; Slater et al., 2019) and conduction velocities of ∼8 m/s (Caminiti et al., 2013; see red cross in Figure 4A). The proposed model yields θ∼0.05 μm for such combination of g_MRI and V, equivalent to a mean axonal radius of 0.45 μm (Figure 4B). This value is consistent with histological analyses of such white matter fibers, with a narrow range of axonal radius (mean∼0.48 μm; Caminiti et al., 2009). For visual transcallosal white matter (MRI-measured g-ratio ∼0.72; conduction velocity ∼10 m/s; Caminiti et al., 2013; Mohammadi et al., 2015; Slater et al., 2019), the proposed model yields θ∼0.23μm (see black cross in Figure 4A), leading to a mean axonal radius of ∼0.63 μm (Figure 4B). This value is also consistent with histological analyses of such white matter fibers, with a broad range of axonal radius (mean∼0.62 μm; Caminiti et al., 2009). For both types of white matter tracts, fibers with a radius above 1.5 μm represent less than 4% of the total number of fibers, consistently with previous histological studies (Aboitiz et al., 1992; Caminiti et al., 2009; Liewald et al., 2014). The value of the model parameter β for frontal and visual transcallosal white matter were 0.68 and 0.73 μm^−α, respectively (Figures 4C,D), in line with histological analyses of the axonal g-ratio in the genu and splenium of the CC (Stikov et al., 2015).

Figure 5 shows the variability of the θ and β estimates in the presence of noise in the in vivo data, computed from combinations of conduction velocity and MRI-measured g-ratio across a plausible range with 700 g_MRI samples and one V sample. The highest errors in θ (∼35%) are found for the smallest values of θ, with only a small effect of the parameter β. The highest errors in β (∼1%) are obtained from a combination of small values of θ and large values of β. Supplementary Figure 3 shows the dependence of the variability of the θ and β estimates on the number of samples of V and g_MRI (β = 0.71 μm^−α and θ = 0.22 μm). With only one estimate of V, as was the case in this study, errors of up to 11% on the θ estimates are observed. The errors on the β estimates are mostly driven by the number of g_MRI samples included in the estimation: for 100 g_MRI samples or more, these errors are below 1% regardless of the number of V samples.

Figure 6 shows the bias in the θ and β estimates arising from an estimation of conduction velocity from a cohort of participants. The bias in θand β is on average 12 and 0.90% across participants, reaching up to 40 and 3%, respectively.

At the sensor level, the group-averaged visual evoked response revealed a positive peak activation, between 115 and 134 ms post-stimulus onset, contralaterally to the stimulus presentation (Figure 7, top). The latency of this positive activation corresponds to the expected latency of the P100 component (Di Russo et al., 2001) and was followed by asymmetrically distributed voltage topographies with maximal voltage values at posterior sites at latencies starting at approximately 153 ms post-stimulus onset (Figure 7, bottom). After approximately 180 ms, post-stimulus onset voltage topographies showed an ipsilateral positivity to the stimulus presentation, suggesting that the P100 activation has traveled to the opposite hemisphere.

Qualitative observations at the electrode-level were confirmed by the source reconstruction results. Group-averaged absolute CD exhibited a sharp increase at approximately 100 ms post-stimulus onset in the right hemisphere compared to the baseline, peaking at 141 ms (Figure 8B for an overview of the spatial distribution of the CD for an exemplar subject). The group-averaged absolute CD on the left visual hemisphere followed the right visual hemisphere, peaking later on at 152 ms.

The IHTT group estimation yielded a value of 11.72 ms with a standard deviation of 2.87 ms. The identified peaks were statistically significantly larger upon visual stimulation when compared to the baseline with p-values of p < 0.01.

Using the estimated group IHTT, we calculated the conduction velocity for the visual transcallosal tract using each participant’s tract length. The average conduction velocity was 13.22 ± 1.18 m/s across participants (Table 1).

Estimates of the model parameters θ and β were computed using the proposed model, from the in vivo samples of the MRI-measured g-ratio along the visual transcallosal tract and estimates of conduction velocity (Table 1). The average value of θ was 0.40 ± 0.07 μm across all subjects and ranged between 0.31 and 0.54 μm (Figure 9A). A θ value of 0.40 μm is equivalent to a mean axonal radius of 0.80 μm, inline with previous estimates from histological studies (0.62 μm; Caminiti et al., 2009). For such a value of θ, axons with a radius above 2 μm represent only < 5% of the total fiber count, in agreement with histological studies that show that axonal radius does not exceed ∼1.5–3 μm in the human brain (Aboitiz et al., 1992; Caminiti et al., 2009; Liewald et al., 2014). The average value of β was 0.67 ± 0.02 μm^−α across all subjects, and ranged between 0.64 and 0.70 μm^−α(Figure 9B). Figure 9C shows the axonal radius distribution P(r) for a representative subject (θ = 0.43 μm), with a confidence interval of 0.27–0.69 μm obtained from the IHTT estimation. The estimated value of β for this subject was 0.68 μm^−α, with a confidence interval between 0.64 and 0.71 μm^−α (Figure 9D).

The proposed model is based on an explicit link between the data acquired in vivo and a limited set of histological properties of white matter axons. The MRI-measured g-ratio is expressed as a function of the axonal radius distribution [P(r)] and the g-ratio of axonal fibers [g(r)], inline with recent studies conducted using MRI and histology data (Stikov et al., 2011, 2015; West et al., 2016). Similarly, axonal conduction velocity — computed from the IHTT estimates — is an ensemble average across the same distribution P(r), assuming an equal contribution of all axons to the EEG data. As a result, both types of in vivo data depend on the same properties of axonal fibers: P(r) and g(r). This approach is supported by recent findings which show that from the numerous histological determinants of conduction velocity (e.g., axonal radius, g-ratio, the conductance of ion channels, diameter and length of Ranvier nodes and internodes), the properties of axons that bring the largest contribution to the determination of conduction velocity are measurable with MRI (Drakesmith et al., 2019).

The mathematical definitions of P(r) and g(r) are grounded on well-established histological findings. P(r) is assumed to follow a gamma distribution, as commonly posited by models of axonal radius distribution (Assaf et al., 2008; Sepehrband et al., 2016). g(r), expressed using a power law, shows a high level of agreement with histological studies (Ikeda and Oka, 2012; Gibson et al., 2014; Figure 2). Besides supporting the proposed model, the histological basis for the expressions of P(r) and g(r) allows freedom in the choice of the model parameters estimated from the data, because they reflect different properties of axon populations that can be set constant or variable according to their relevance in a neuroscience application of this model.

From the set of 4 model parameters, we opted to set the mode M of the axonal radius distribution to a constant value, based on the histological literature (Tomasi et al., 2012; Liewald et al., 2014). Similarly, the parameter α was set from a calibration with histological findings (g_REF) (Ikeda and Oka, 2012; see section “Axon Morphological Properties”). The estimated parameters — θ, the right tail of P(r) and β, the scaling parameter of g(r) — enabled the simultaneous assessment of axonal radius and myelination, across the distribution of axons in the visual transcallosal white matter tract. The average value of θ was 0.40 μm across participants, leading to a mean axonal radius of 0.80 μm, inline with previous estimates from histological studies (0.62 μm; Caminiti et al., 2009). The average value of β was 0.67 μm^−α across participants, consistent with histological measures of the axonal g-ratio in the splenium of the CC (Jung et al., 2018).

The proposed approach inherits the limitations of the MRI methodologies used to estimate the intra-cellular and myelin volume fractions, and subsequently the MRI-measured g-ratio. In the current application, the NODDI model was used to estimate the intra-cellular volume fraction (Zhang et al., 2012). This model assumes identical parallel diffusivity in the extra- and intra-cellular spaces, set to 1.7 μm²/ms (Zhang et al., 2012). This simplifying assumption might lead to potential bias of the parameter estimates (Jelescu et al., 2015). We estimated that a change in diffusivity within a realistic range (1.5–1.9 μm²/ms, Guerrero et al., 2019) leads to a change of 4–5% and 1–2% for θ and β, respectively, smaller than the bias arising from the use of group-averaged conduction velocities (12% for θ and 0.90% for β). Alternatives models (e.g., Assaf and Basser, 2005; Fieremans et al., 2011; Campbell et al., 2018; Ellerbrock and Mohammadi, 2018) may be considered in light of their assumptions, as well as their applicability given the available data.

The IHTT was estimated from the group averaged CDs, which allowed for easy identification of the first two maxima of activation in the two hemispheres (Figure 8). The IHTT derived as a result of the time interval separating the peaks of CDs at the group level (11.72 ± 2.87 ms) falls within the range (i.e., 8 and 30 ms) of previous IHTT estimates based on a priori selection of voltage measurement at electrodes at occipital sites (e.g., Saron and Foxe, 2003; Westerhausen et al., 2006; Whitford et al., 2011; Friedrich et al., 2017; Chaumillon et al., 2018). In our study, we opted for a CD-based estimation of the IHTT as the closest reflection of the evoked neural activity in the regions of interest, consistent with the corresponding white matter tract selection (section “MRI-Based Estimation of the g-ratio”). In addition, this approach helps to overcome the ambiguity of electrode selection in electrode-based IHTT estimations. The resulting estimates of conduction velocity are in agreement with the values reported by Caminiti and colleagues (10 m/s, Caminiti et al., 2013), obtained from human histological samples of the same fiber tract. Our results are also consistent with a more recent report that uses a different approach based on Bayesian networks to map the flow of information following left visual stimulation (Deslauriers-Gauthier et al., 2019). The authors observed a transfer of information from the right to the left occipital cortex, between 140 and 160 ms (Deslauriers-Gauthier et al., 2019), in agreement with the latencies of the right and left activations observed in our analysis.

Ideally, IHTT estimates should be performed at the single-subject level to provide microstructure measures specific to each subject. However, the presence of multiple peaks of activation in some individual datasets and the resulting ambiguity in defining the maxima of activation consistently across subjects prevented the estimation of individual IHTT values. These singularities might arise from differences in brain morphology (Saron and Foxe, 2003), inaccurate estimation of the sources, or from the lack of an objective criterion for selecting subjects fulfilling our assumptions on the expected pattern of activations. The difficulty in estimating the IHTT at the subject level has been pointed out in many other studies that showed inconsistent IHTT values across participants and counterintuitively, even negative values in some cases (e.g., Saron and Davidson, 1989; Marzi et al., 1991; Westerhausen et al., 2006; Friedrich et al., 2017). Estimation of the IHTT might be improved with further development on the computation of the inverse solution (Plomp et al., 2010; Mahjoory et al., 2017), and by introducing priors to constrain the time courses of the activations.

The proposed model is based on the combination of a structural measure of white matter (MRI-measured g-ratio) and a measure of brain function (axonal conduction velocity). The validity of this model relies on the assumption that both measures may be obtained for a given white matter tract. Anatomical delineation of a white matter tract is generally a routine procedure thanks to MRI tractography techniques (Caminiti et al., 2013; Horowitz et al., 2015; Tournier et al., 2019), allowing the sampling of the MRI-measured g-ratio data along this tract (Schiavi et al., 2022). On the other hand, the underlying mechanisms and anatomy of the inter-hemispheric transfer of the evoked visual activity with Poffenberger’s paradigm are still under investigation. Previous literature has used this paradigm to demonstrate an increase in IHTT in acallosal subjects, highlighting the primary role of the CC in visual interhemispheric transfer (Marzi et al., 1991; Di Stefano et al., 1992; Aglioti et al., 1993; Tassinari et al., 1994). In addition, Westerhausen et al. (2006) showed that IHTT is significantly correlated with the structural integrity of the posterior CC, suggesting splenium fibers as the most likely pathway for visual interhemispheric transfer.

The source and target cortical areas of the visual interhemispheric transfer remain to be fully identified. Previous studies have demonstrated the existence of a small patch of transcallosal axons between visual areas 17 and 18 (Clarke and Miklossy, 1990; Aboitiz and Montiel, 2003). This motivated our choice of source and target cortical areas, which led to conduction velocity estimates consistent with previous studies (Caminiti et al., 2013). However, it has been suggested that these connections alone might be insufficient to produce an effective interhemispheric transfer (Innocenti et al., 2015) and higher visual processing areas might be involved.

The proposed model requires a measure of the MRI-measured g-ratio and axonal conduction velocity in a specific white matter tract of interest. Since measurements of conduction velocity may only be conducted for a limited set of white matter tracts in the human brain, this model, in its current form, cannot be extended to the entire white matter in contrast to other models (e.g., Assaf et al., 2008; Zhang et al., 2011). Instead, this model is geared toward the detailed characterization of a restricted set of white matter tracts. In its proposed form, two parameters relating to the morphology of white matter axons may be estimated from a total of four. Extending the number of estimated parameters requires the use of additional data integrated into the proposed framework. The nature of this data needs to be carefully considered. Alternatively, as was outlined in this work, a subset of the model parameters may be set to a constant value based on histological studies. A prime candidate is the mode of the axonal radius distribution (M), preserved between white matter tracts, individuals, and animal species (Tomasi et al., 2012; Liewald et al., 2014). This choice may need to be reconsidered for the study of brain pathologies that might differentially affect axons of different sizes. Other potential candidates are the parameters α or β that describe the dependence of the fiber g-ratio on axonal radius. In this study, we opted to set α to a constant in expectation of fiber myelination differences identical for all radii. Alternatively, setting β to a constant value may be preferred to allow for the estimation of α, in cases where myelin thickness differences are expected to depend on axonal radius.

In the current study, the estimate of the IHTT was calculated from the average of the EEG-based CD across participants. This led to an average bias in θand β of 12% and 0.90% across participants, reaching up to 40 and 3%, respectively (Figure 6). While small compared to differences between tracts (∼350% for θ and ∼7% for β, see Figure 4), this is of the order of the inter-subject differences for these parameters (Figure 9) and prevents the estimation of axonal morphological features at the individual level. This limitation represents the primary avenue for future improvements. The question of the accuracy of the parameter estimates and their comparison with histological data might arise subsequently. Estimation of the IHTT from alternative techniques such as transcranial magnetic stimulation might also be considered (Lo and Fook-Chong, 2004; Spitzer et al., 2004; Basso et al., 2006; Deftereos et al., 2008; Marzi et al., 2009). We highlight that the parameters of the proposed model touch on properties of brain tissue that have received little attention in histological studies to date. These include the radius dependence of the axonal g-ratio and the radius dependence of fiber myelination change in health and disease. Validation of the proposed model may therefore bring opportunities for new research avenues for histological studies of the human brain.