hMRI – A toolbox for quantitative MRI in neuroscience and clinical research

1Quantitative MRI (qMRI) finds increasing interest in neuroscience and clinical research because it is not only more sensitive, but also more specific, to microstructural properties of brain tissue such as axon, myelin, iron and water concentration than conventional weighted MRI. In conventional weighted MRI, the image grayscale values have arbitrary units and the value in a given voxel will depend on a large number of factors, such as the sequence type (e.g. the magnetisation-prepared rapid gradient echo, MPRAGE versus modified driven equilibrium Fourier transform, MDEFT for T1-weighted anatomical images), sequence parameters (e.g. repetition time,, echo time,, or flip angle), and hardware effects (e.g. transmit and receive profiles and any scaling factors). In addition, the value will depend on multiple physical tissue properties such as the longitudinal and transverse relaxation times, and, or the proton density,. qMRI accounts for these varied effects in order to increase the specificity of the estimated metrics and eventually quantify specific physical tissue properties. In qMRI, the estimated physical value has a direct meaning and is quantified in standardised units (e.g. in seconds). This standardised nature further increases the comparability across sites and time points, which may improve the sensitivity of multi-site studies and longitudinal analyses of development, plasticity and disease progression. A biophysical interpretation of physical qMRI parameters or a combination of qMRI with biophysical modelling (e.g.,, or ) enables the in vivo characterisation of key microscopic brain tissue parameters, which previously could only be achieved with ex vivo histology. This concept is called in vivo histology using MRI (hMRI, ).

2The estimation of quantitative and semi-quantitative metrics commonly includes one or more of the effective transverse relaxation rate, the longitudinal relaxation rate, the proton density, the magnetisation transfer saturation and a number of diffusion MRI (dMRI) metrics. However, the majority of fundamental and clinical neuroscience studies either refrain from acquiring qMRI data or the quantitative approach is limited to dMRI only. One reason for this might be that standardised qMRI imaging protocols (such as the protocol for the Human Connectome Project) and processing software (see a summary in ) are readily available for dMRI but less so for other qMRI techniques.

3Consequently, the neuroscience and clinical research community lacks a standardised qMRI implementation to handle the wide diversity of data acquisition types and estimate parameters such as,, and saturation, which are sensitive to iron, myelin, and water content in tissue microstructure, and thus provide complementary information to the axonal properties revealed by dMRI. For example, is often estimated using gradient recalled echo data acquired with multiple echo times to create high resolution maps that show strong contrast between different types of brain tissue. Similarly, and can be derived from two acquisitions varying the excitation flip angle. For accurate estimation of the latter qMRI parameters, an implementation must adequately correct for instrumental biases such as inhomogeneous transmit and receive fields. Such correction should be based on the additional mapping of these fields, but the implementation should also provide solutions when these fields have not been measured. For example in clinical settings, MR sequences for estimating instrumental biases are often unavailable due to time, hardware or software constraints. For such studies, the facility to correct for instrumental biases retrospectively using image processing methods that do not rely on additional MRI acquisitions is highly desirable.

4A number of open-source tools have been developed to support qMRI use and make qMRI more broadly accessible to neuroscience and clinical research (see for example qMRLab (, https://github.com/qMRLab/qMRLab), QUIT (, https://github.com/spinicist/QUIT), mrQ (, https://github.com/mezera/mrQ) and QMAP (https://www.medphysics.wisc.edu/∼samsonov/qmap/). These tools include, to various extents, data acquisition guidelines, tools for protocol simulation and optimisation, a multitude of models and methods for data fitting and estimation of qMRI parameters, and visualisation tools.

5While most of these tools focus on model fitting and generation of qMRI maps, the question of spatial and statistical processing of these qMRI maps is not directly addressed. For spatial processing and group level statistical analysis, the established tools are primarily designed for diffusion MRI (see for example TBSS in FSL and TRACULA in FreeSurfer). In addition, a number of custom made tools have been developed based on established neuroimaging software, e.g.: applications of the FreeSurfer surface projection software to compare quantitative relaxation and susceptibility data on the cortex, or usage of the VBM framework to process qMRI data across the whole brain on a voxel-by-voxel basis ( and for dMRI metrics, and Table 1 for,, and saturation). One challenge common to all these tools is to find a proper method to locally preserve the qMRI parameters after (non-linear) spatial registration. The methods using the statistical parametric mapping (SPM) framework typically reduce residual misalignment between images by isotropic spatial smoothing. However, applying this framework directly to qMRI data would introduce partial volume effects at tissue boundaries and corrupt the quantitative values. Hence a modified smoothing approach, which aims to achieve within class smoothing only, is preferred. The majority of the studies in Table 1 took advantage of voxel-based quantification (VBQ), an approach introduced by and developed for the comprehensive multi-parameter mapping (MPM) approach to correct for potential error introduced by spatial smoothing.

6A particular instantiation of qMRI developed at 3T, the MPM approach, spans data acquisition, modelling and bias correction of three multi-echo spoiled gradient echo volumes to generate,,, as well as semi-quantitative saturation maps. This framework enables time-efficient whole brain mapping with high isotropic resolution of 800 μm in 27 min or reduced MPM protocol (no MT saturation) at 1 mm isotropic resolution in 14 min at 3T and has even enabled the acquisition of ultra-high-resolution quantitative maps with 400 μm resolution at 7T. This framework has been used in a variety of fundamental and clinical neuroscience studies (Table 1) focussing on: (a) improving the segmentation of deep grey matter structures, (b) evaluating the myelin and iron concentration in the brain and spinal cord, (c) linking structure and function in the cortex, and (d) using the MPM parameter maps as proxies for biophysical tissue models.

7In this paper, we present the hMRI-toolbox, a comprehensive open-source toolbox that streamlines all the processing steps required to generate,, and saturation maps and provides appropriate spatial processing for group analyses. The flexible nature of the toolbox makes it applicable to a wide range of data types, from the full MPM protocol to subsets of it, including single contrast echo trains for mapping or variable flip angle data for mapping of and using multi-echo or even single-echo data. The toolbox is embedded in the SPM framework (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/), profiting from the highly accurate spatial registration into a common space and the variety of established statistical inference schemes. The spatial processing part of the toolbox can be applied to any set of rotationally-invariant qMRI maps, including a number of diffusion MRI parameters and all common qMRI metrics.

8The MPM multi-echo protocol was introduced in and for estimating the longitudinal relaxation rate, the effective transverse relaxation rate, the proton density and the magnetisation transfer and generalises a number of acquisition protocols. It typically involves acquiring six to eight images at different echo times for each of the PD-, T1-and MT-weighted acquisitions in an RF and gradient spoiled gradient echo sequence (referred to as T1w, PDw and MTw echoes, respectively). The hMRI-toolbox can flexibly deal with a large range of site- and study-specific acquisition schemes, from the full MPM protocol to subsets of it, including single contrast echo trains for mapping or variable flip angle data for mapping of and using multi-echo or single-echo data, comparable to e.g. DESPOT1. Example MPM acquisition protocols can be found at http://hmri.info.

9We give here, and in Fig. 1, a short overview of the theory underlying the qMRI map creation process. For a more detailed outline of the theory and the estimation procedures applied in the hMRI-toolbox see Appendix A.

10The signal from the multi-echo PDw, T1w and MTw acquisitions can be described by the Ernst equation. The effective transverse relaxation rate can then be derived from the dependence of the signal. The unified description of the multi-echo data from all three contrasts into a single model, denoted as ESTATICS, provides a more robust estimation of with a higher signal-to-noise ratio compared to separate estimations (Fig. 1a). Using approximations of the signal equations for small repetition time and small flip angles α, the longitudinal relaxation rate, the apparent signal amplitude map (proportional to the proton density ) and the magnetisation transfer can be estimated. At this point (Fig. 1b), the generated maps are biased by transmit (Fig. 1c) and receive (Fig. 1d) field inhomogeneities. The hMRI-toolbox provides correction methods for these bias fields based on specific transmit and receive field measurements or image processing methods. While influences the local flip angle and hence all three maps are affected, the RF sensitivity bias field only influences the map (in the absence of subject motion).

11The toolbox can also handle the situation where only a subset of data is available. For example,, and can still be estimated when no MTw acquisitions are acquired, alone can be estimated when neither MTw nor T1w acquisitions are available (single multi-echo PDw data)., and saturation maps can be generated from single echo PDw, T1w and MTw images, not requiring multi-echo acquisitions. The theory and map creation tools also encompass the creation of maps from other variable flip angle approaches, such as DESPOT1 or maps from multi-echo data, such as certain susceptibility mapping/weighted imaging approaches.

12The latest version of the toolbox can be downloaded from the hMRI-toolbox page (http://hmri.info) as a zip file (containing the last official release) or by cloning the git repository (https://github.com/hMRI-group/hMRI-toolbox) to keep up-to-date with the latest incremental developments.

13Updated documentation is available as a Wiki (https://github.com/hMRI-group/hMRI-toolbox/wiki). It includes installation instructions, an example dataset, a tutorial and a detailed description of the implemented functionalities. Information on releases and versioning, development and contribution guidelines are also provided.

14The toolbox has been developed and tested with MATLAB versions 8.0 (R2012b) to 9.3 (R2017b) and SPM12 from version r6906 onwards. Since the hMRI developments will be synchronised with SPM developments, it is recommended to use the most recent SPM release (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) to benefit from the latest developments.

15The hMRI-toolbox is free but copyright software, distributed under the terms of the GNU General Public License as published by the Free Software Foundation (as given in file Copyright.txt). Further details on “copyleft” can be found at http://www.gnu.org/copyleft/.

16In particular, the hMRI-toolbox is supplied as is. No formal support or maintenance is provided or implied. Since the toolbox was developed for academic research, it comes with no warranty and is not intended for clinical use.

17An MPM example dataset from a healthy subject for demonstrating the hMRI-toolbox features was acquired on a 3T Prisma system (Siemens Healthcare, Erlangen, Germany) at the Wellcome Centre for Human Neuroimaging, London, UK. This dataset can be downloaded from http://hmri.info and is described in detail in.

18The hMRI-toolbox is organised into five main modules (Fig. 2): Configure toolbox, DICOM import, Auto-reorient, Create hMRI maps and Process hMRI maps. While the Configure Toolbox and Process hMRI maps modules can be run for a group of subjects, the DICOM Import, Auto-reorient and Create hMRI maps modules must be run for each subject and each session separately (if several datasets acquired per subject, e.g. in a longitudinal study). A brief description of each module is provided below and details can be found in the Appendices.

19The hMRI-toolbox provides the user with a set of default acquisition and processing parameters for most common acquisition protocols without any further customisation. However, customisation is possible and necessary to broaden the toolbox usability to a wider range of protocols, scanners and vendors. Therefore, the Configure toolbox module allows the user to select either standard or customised default parameters, to match the user's own site- or protocol-specific setup and to be used in the subsequent modules (see details in Appendix C.1).

20DICOM import is a tool to convert DICOM data into NIfTI files. During conversion, the whole DICOM header is stored as JSON-encoded metadata in a file along side the NIfTI images. This feature is implemented in SPM12 from release r7219 (November 2017), following the hMRI-toolbox implementation. Note that Philips-specific rescaling factor is applied at conversion, starting from version v0.2.0 of the hMRI-toolbox and release r7487 of SPM12. The hMRI-toolbox also provides metadata handling functionalities to retrieve parameter values and store processing parameters in the JSON-encoded metadata file (see Table 3 for an example of processing information stored as metadata). Detailed description of the DICOM import module, the metadata handling tools and BIDS compliance aspects is provided in Appendix B.

21The reorientation of the images towards a standard pose, setting the anterior commissure at the origin and both anterior and posterior commissure (AC/PC) in the same axial plane, as defined in MNI space, is a common step that increases the consistency in individual head positions prior to normalisation or segmentation. For example, SPM's segmentation is sensitive to the initial orientation of the images. Therefore, Auto-reorient provides a simple tool for automatically and uniformly reorienting a set of images prior to any further processing including multi-parameter map calculation. Note, that the reorientation modifies the orientation information in each image header, but the reoriented images are not resliced. For more details, see Appendix C.2.

22The Create hMRI maps module computes quantitative as well as semi-quantitative estimates of,, and saturation from unprocessed multi-echo T1w, PDw and MTw spoiled gradient echo acquisitions. The map creation module corrects the qMRI estimates for spatial receive (Section 3.7.3) and transmit (Section 3.7.2) field inhomogeneities.

23The module takes the (possibly reoriented) series of multi-echo spoiled gradient echo images as input for the creation of quantitative as well as semi-quantitative maps (Fig. 3) of,, and saturation as described in the Background Section 2 and Appendix A. The number and quality of the output maps depends on the contrasts (PDw, T1w, MTw) and number of echoes available, and on the availability of additional bias field measurements. A single multi-echo PDw contrast allows for the calculation of a single map. If T1w images are additionally provided, and maps will also be generated. The map estimation requires the acquisition of additional MTw images. The map creation also involves optional correction for ( transmit) bias field (Section 3.7.2) and the (RF receive sensitivity) bias field (Section 3.7.3) as well as spoiling imperfections.

24The map creation module includes the determination of transmit bias field maps ( expressed in p.u. Of the nominal flip angle) for transmit bias correction of the quantitative data. Several methods are implemented. Depending on the choice of the specific method the GUI requires the user to provide adequate input files. Further details on the supported correction methods can be found in Appendix C.3 and the respective original publications.

25Three options are available to correct for RF receive sensitivity bias within the Create hMRI maps module. Two of them rely on measured RF sensitivity maps (Single or Per contrast options) while the third method is data driven (Unified Segmentation option: no input sensitivity map required). Although not recommended, it is possible to disable RF sensitivity correction altogether by selecting the None option. While options Single and Unified Segmentation assume that the sensitivity profile is consistent between contrasts (i.e. small inter-contrast subject movement is assumed), the Per contrast option accounts for inter-contrast variation in RF sensitivity profile due to larger subject motion. Details on the different RF sensitivity bias correction methods can be found in Appendix C.4 and in the respective publications.

26By default, the estimated quantitative maps are output into a Results subdirectory within the folder of the first PDw echo. Alternatively, a user-defined folder for the output of the toolbox can be selected in which the Results directory will be created. The estimated qMRI maps are saved in the Results directory, with supplementary files are output in the Results/Supplementary subfolder. The basename for all qMRI maps is derived from the first echo of the PDw image series, see Table 2 for brief description. If data is reprocessed, a new sub-folder is created.

27The Process hMRI maps module provides dedicated tools and tissue probability maps for the spatial processing of quantitative MRI maps based on the corresponding SPM framework. The spatial processing pipeline for hMRI data relies on three main operational steps (Fig. 4): (1) segmentation, (2) non-linear spatial registration into common space and (3) tissue-weighted smoothing, using three different sub-modules that are further detailed in the Appendix, Table 4. Furthermore, a fully integrated processing pipeline is provided as an additional sub-module to facilitate standard data processing without the need to combine the individual steps in this module. Details on the three sub-modules and the integrated pipeline are provided in Appendix C.5.

28The standard SPM statistical analysis and modelling approaches such as mass-univariate General Linear Modelling can be applied to the spatially processed maps, see, e.g., and. Additionally, the multiple parameter maps lend themselves to multi-variate analyses approaches as well.

29This paper introduced the hMRI-toolbox, which is embedded in the SPM framework and allows for the estimation and processing of four quantitative MRI parameter maps: the longitudinal and effective transverse relaxation rates and, the proton density and the (semi-quantitative) magnetisation transfer saturation. This introduction includes a comprehensive summary of the MPM signal model as well as the currently available correction methods for the various bias sources that, if not corrected for, might impair the quantification. Finally, the processing steps for a dedicated SPM analysis of quantitative MRI maps (denoted as the VBQ approach) were presented, correcting for the potential partial volume effects introduced by spatial smoothing.

30The name of the toolbox (h-MRI) originates from the concept of in vivo histology of tissue microstructure using MRI. Hereby, the quantitative parameter maps generated with this toolbox provide key input parameters for biophysical models that are designed to non-invasively estimate specific microstructural tissue properties (see Table 1 for example studies).

31Considerations and interpretation of the MPM approach. While the signal model (see Appendix A for details) used in this toolbox is based on the Ernst equation and thus provides a comprehensive means of calculating a set of physical quantitative (and semi-quantitative) parameters, including,, and, we would like to emphasise that more sophisticated models can be derived to relate the parameters more directly to the underlying biophysical mechanisms and tissue characteristics.

32The contrast provided by the metric is associated with different realisation of iron deposits (for a review, see e.g. ), myelination of axons, their orientation and chemical exchange. Multi-compartment models can be used, e.g., to separately describe the orientation-dependence of myelinated fibre pathways in parameter maps (e.g.,;; ).

33At a given field strength, the contrast is determined by the micro-structural tissue properties such as the local mobility of water molecules, the macromolecular content and the local concentration of paramagnetic ions such as iron or gadolinium-based contrast agents. It has been shown that depends within limits linearly on these tissue properties. The and metrics can be used to describe the dependence of on these components.

34saturation is a proxy measure of the bound-pool water fraction. It provides information about the macromolecular content of the micro-structural environment and is often used as a marker for myelin content, see e.g. and. Under equivalent conditions for the off-resonance pre-pulse the same saturation values are expected. However, this will not be the case if the properties of the pre-pulse are changed across measurements. Therefore, we refer to the saturation measure as being semi-quantitative. The saturation map differs from the commonly used ratio (MTR; percent reduction in steady state signal) by explicitly removing the bias introduced by the spatially varying relaxation time and -transmit field. Additional minor corrections for transmit field inhomogeneity in the maps were applied as described in. The reduced spatially varying bias leads, e.g., to a higher contrast in deep brain structures than MTR and to reduced variance in the data. Note that the saturation measure does not only depend on the bound-pool fraction but also on the exchange between the bound and free water pools (see e.g. ). A more direct measure of the bound-pool fraction is provided by quantitative, which requires more time-consuming data acquisition, typically limiting the map's spatial resolution (e.g., 2 mm isotropic in ).

35Sources of bias in qMRI and limitations to bias corrections. Correction of instrumental characteristics and artefacts is an essential prerequisite for quantitative MRI. Sources of biases and artefacts include primarily transmit and receive fields, imperfect spoiling, -bias and head motion. The artefact correction methods provided in the hMRI-toolbox are highly flexible, offering solutions to process reference measurements (e.g. transmit/receive field measurements carried out with a number of customised or product sequences) to correct for instrumental artefacts, as well as achieve optimal results even when no adequate measurements are available.

36The ideal imaging protocol includes dedicated measurements of transmit and receive field inhomogeneities, which are typically based on customised sequences (Fig. 5a–d for RF receive bias correction, and Fig. 6a for transmit bias correction). When customised sequences are not available, a standard spoiled gradient echo product sequence can be used to acquire data for low-resolution receive field mapping. Similarly, when not available, customised reference and calibration sequences such as 3D_EPI or 3D_AFI for correction of transmit effects can be replaced by manufacturer's service sequences such as TFL_B1_map or RF_map (as examples on Siemens scanners). The toolbox provides the option to process data from several different transmit bias field mapping techniques or to use transmit bias field maps pre-computed outside the toolbox (see section 3.7.2 and Appendix C.3).

37The measured transmit and receive fields can be affected by diverse sources of error leading to imperfect corrections. For example, residual misalignment between measured receive (or transmit) field and the spoiled gradient echo images can be one reason for such imperfections. In particular, when between-contrast (PDw, T1w, MTw) motion is large, discrepancies between head position for a single receive (or, to a lesser extent, transmit) field measurement and head position for all or some of the spoiled gradient echo images lead to additional motion-related bias in the quantitative maps (Fig. 5). In such a case, a per-contrast RF receive sensitivity measurement is preferable and can account for the between-contrast dynamic variation (see Appendix C.4, Receive field sensitivity measurements and Fig. 5). However, measured RF receive sensitivity maps such as described in Appendix C.4 can also suffer from residual modulations by the receive field of the body coil, which serves as a reference and whose inhomogeneity is not accounted for. Such modulation cannot be directly corrected for using the measured transmit field of the body coil at 3T and higher fields due to the non trivial applicability of the reciprocity principle at such field strengths. As a result, when no large between-contrast motion is observed, RF sensitivity bias correction using the data driven receive field estimation (described in Appendix C.4) may prove more effective altogether. Such a data driven method could also be applied (with specific optimisation of the US regularisation parameters to the body coil's receive field profile) to correct for the above residual body coil receive field modulation. Finally, at 7T where RF body coils are not available, the currently implemented RF sensitivity measurements and bias correction approach are not applicable. All the above aspects are active fields of investigation, optimisation and validation that are also an integral part of the future developments of the toolbox. Similarly, transmit field mapping techniques can be inaccurate. For a comparison of the frequently used -transmit field mapping techniques and a description of their respective accuracy and sources of uncertainty we refer to and.

38When no transmit and/or receive field inhomogeneity maps have been measured, which often happens in clinical settings due to time constraints, the toolbox provides the option to use image processing methods based on the Unified Segmentation approach for transmit bias correction (UNICORT,; see Fig. 6b) or RF sensitivity bias correction (Fig. 5e). The Unified Segmentation approach takes advantage of the fact that bias corrections can be applied post hoc in good approximation for small read-out flip angles and short. This requires no additional acquisition time but produces quantitative maps of lesser accuracy with some residual receive and/or transmit field modulation compared to a correction with measured references. Note that the Unified Segmentation approach, whether applied for transmit (UNICORT) or RF sensitivity bias correction, has been optimised for the Siemens TIM-TRIO MR system using the body RF coil for transmission and the 32-channel receive head coil. The corrections will perform appropriately for coils with similar transmit or receive field profiles, but might require further adjustments otherwise (see for UNICORT optimisation). For more details, see and the hMRI wiki (the latter also provides further information on customized usage).

39The proposed MPM protocol uses RF and gradient spoiling to minimise undesired transverse net magnetisation. Imperfect spoiling, which depends on the precise sequence protocol settings, can leave a residual bias in the map if no further correction is used. For specific MPM protocols using the customised sequences, the hMRI-toolbox provides a correction for imperfect spoiling, see Eq. (A.15) in Appendix A.5. By default this correction is disabled but can be enabled through the toolbox customisation provided by the Configure toolbox module (Appendix C.1 and the hMRI-toolbox Wiki).

40The estimation of maps can be biased by relaxation effects if not accounted for (Eq. (A.4) in Appendix A). Two correction methods, based on extrapolation of the data to or relying on the estimated maps respectively, are implemented in the toolbox. These correction methods require a number of echoes to be acquired for a robust fit of the exponential decay to derive magnitude images and maps. In the case of a single-echo variable flip angle dataset, the -weighting correction cannot be applied and the estimated and (to a much lesser extent) the and maps will be biased by -modulations. Also, maps are biased in areas with severe susceptibility artifacts. Note that the JSON metadata file associated with the respective PD parameter map contains information about the processing steps and thus of potential modulation.

41Head motion is widely recognised as a major source of artefacts in MR images, with severe consequences for quantitative MRI and morphological measures of the brain. While quantitative measures of image quality have been introduced, visual inspection remains the most common means of rating data quality despite its limited sensitivity and inter-rater variability. The hMRI-toolbox provides summary measures of head motion within and between the acquisitions of each image volume (intra- and inter-scan motion) ( and Appendix D). The provided index of intra-scan motion has been tested against the history of head motion, recorded in real-time during the scans. Note that these intra- and inter-scan head motion measures could potentially be combined to guide toolbox users to objectively classify their data according to quality, for example to exclude or downweight poor-quality data of individuals in a statistical group analysis.

42Spatial processing pitfalls. Since spatial processing in the hMRI-toolbox is embedded in the SPM framework, it is subject to the same limitations as any typical VBM study, including spatial normalisation accuracy, segmentation errors and partial volume effects. Auto-reorient is an option that can help improve the segmentation. However, it has to be done carefully. Poor signal-to-noise ratio, contrast-to-noise ratio, or outliers in the MPM input images may impair the reorientation procedure. In general, it is good practice to visually inspect the results of the hMRI pipeline to detect any obviously suspicious results. Piloting the processing pipeline using a batch for a single healthy subject dataset, combined with a careful check of the log files (in Results/Supplementary, metadata including summary description as shown in Table 3), and comparison with the computed MPM maps from the sample dataset, is recommended.

43Residual misalignments between qMRI maps of individual participants in common space will be present despite the high degree of spatial correspondence that can be achieved by the non-linear warping algorithms available in SPM, e.g., by DARTEL. These misalignments are typically reduced in the VBM-framework by spatial smoothing. To correct for the partial volume effects at tissue boundaries that can be introduced by spatial smoothing, the hMRI-toolbox provides the dedicated smoothing approach that has been described in Fig. 4.

44The choice of the appropriate smoothing kernel and its performance compared to alternative methods (e.g. TBSS or TSPOON) is still a subject of active research. An alternative method for reducing spatial misregistration in the cortex might be surface-based registration algorithms. However, a comparison study by between volume-based and surface-based registration methods could not demonstrate clear superiority of one or the other approach. In fact this is an active area of research to better understand the relative benefits and pitfalls of each approach.

45hMRI-toolbox for different MRI scanner platforms. The first task for interested users of the hMRI-toolbox is the setup of the MPM acquisition protocol. To facilitate standardisation, a set of example protocols for the customised MPM sequences on Siemens platforms as well as standard sequences available on Siemens and Philips platforms are provided on the http://hmri.info website. Those protocols take advantage from the fact that MPM sequences primarily rely on multi-echo spoiled gradient echo sequences that are available on all modern MRI scanners. Even though the experience with implementing the specific MPM protocol on the Philips platforms is limited, first important steps have been achieved, and the MPM framework together with the hMRI-toolbox will be used in a multi-site clinical trial (NISCI trial,,) including Philips and Siemens MR systems.

46Multi-scanner and multi-vendor data sets will require adjustments in terms of data handling and processing with the hMRI-toolbox. Most MPM studies up to now were carried out on Siemens MRI scanners using customised MPM sequences. Consequently, the toolbox is optimised for this scenario. New issues might arise when implementing the MPM protocol and using data from other vendors or conventional product sequences. For instance, different MT pre-pulse implementations will lead to changes of the saturation map, which will require appropriate inter-scanner calibration. Moreover, not every DICOM to NIfTI conversion software appropriately handles image intensity scaling, as reported e.g. for Philips data, leading to spurious intensity differences affecting the quantification. Making the hMRI-toolbox data formats fully BIDS compliant by defining the ontology of acquisition parameters necessary for the creation of the quantitative maps (see Appendix B.5) further supports the use on multiple platforms and vendors, and is a high priority of the ongoing developments.

47Applicability of the toolbox to different MRI platforms also involves ultra-high field MR systems. With the fast progress of ultra-high field scanners (7T and higher) on the MRI market, high-resolution data will become more routinely available and provide access to e.g. laminar-specific information (see 400 μm new generation MPM, ), while also posing new challenges (e.g. field inhomogeneities) that the hMRI-toolbox will have to address.

48Future directions. The hMRI-toolbox has been developed as a scientific collaborative project. As such its developments aim at making it broadly available, capitalising on its flexible and open-source implementation, and adjusting to data sets acquired on multiple MRI platforms (see hMRI-toolbox for different MRI scanner platforms above).

49Sensitivity to small inter-individual changes of microstructure (e.g. plasticity) and the variation of change across subjects (e.g. in development, see https://www.biorxiv.org/content/early/2018/07/26/328146) is another challenge for longitudinal qMRI. To that end, the bias of the qMRI estimates in the presence of motion has been investigated. Retrospective robust estimation of parameters (outlier rejection) has been suggested in and could be implemented in future releases of the hMRI toolbox. Robustness must also be considered in parallel to the spatial resolution versus sufficient SNR level tradeoff to improve the sensitivity of the technique to small developmental changes and plasticity. Thus, spatially adaptive noise removal methods along with appropriate handling of the Rician bias problem are important future improvements of the hMRI-toolbox.

50Currently, quality assessment (QA) is provided in the hMRI-toolbox via a set of indicative parameters that should be used at the user's discretion (Appendix D). Future work will focus on further validating and implementing an automated QA of the raw data and generated maps, providing figures from representative populations and protocols.

51Extensions of the hMRI-toolbox can easily fit within its modular implementation. As short term future additions, the following three modules and extensions are considered. An additional module that efficiently calculates the protocol-specific correction parameters required to account for imperfect spoiling (Appendix A.5) is planned. Quantitative susceptibility mapping (QSM), taking advantage of the existing phase images acquired with the MPM protocol, is a second extension. Finally, as suggested by the ”h” in hMRI, adding new biophysical models that take advantage of the multi-contrast MRI data and generated qMRI maps for in vivo histology is another priority of future developments. An example for such a direct extension of the hMRI-toolbox could be the MR g-ratio model. The MR g-ratio (the ratio between inner and outer diameter of a myelinated axon) is a geometrical microstructural tissue property that can be derived by combining myelin-sensitive qMRI maps from the hMRI-toolbox (e.g. or maps) with the axonal-sensitive maps obtained with existing SPM tools (e.g. the ACID toolbox and the DTI & Fiber tools).

52The hMRI-toolbox is a highly flexible software package that provides a computationally time-efficient, robust and simple framework to generate and process qMRI parameter maps sensitive to myelin, iron, and water content. It profits from the powerful and easy-to-use spatial and statistical analysis tools in SPM, and can be readily combined with existing SPM tools for quantitative estimation of parameter maps sensitive to complementary information such as axonal properties. The ongoing developments address the use of open-science data formats and extensions into biophysical models for direct microstructure mapping. As such, the hMRI-toolbox is a comprehensive and readily extendable tool for estimating and processing qMRI data for neuroscience and clinical research.