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

Neuroscience and clinical researchers are increasingly interested in quantitative magnetic resonance imaging (qMRI) due to its sensitivity to micro-structural properties of brain tissue such as axon, myelin, iron and water concentration. We introduce the hMRI-toolbox, an open-source, easy-to-use tool available on GitHub, for qMRI data handling and processing, presented together with a tutorial and example dataset. This toolbox allows the estimation of high-quality multi-parameter qMRI maps (longitudinal and effective transverse relaxation rates R1 and R2⋆, proton density PD and magnetisation transfer MT saturation) that can be used for quantitative parameter analysis and accurate delineation of subcortical brain structures. The qMRI maps generated by the toolbox are key input parameters for biophysical models designed to estimate tissue microstructure properties such as the MR g-ratio and to derive standard and novel MRI biomarkers. Thus, the current version of the toolbox is a first step towards in vivo histology using MRI (hMRI) and is being extended further in this direction. Embedded in the Statistical Parametric Mapping (SPM) framework, it benefits from the extensive range of established SPM tools for high-accuracy spatial registration and statistical inferences and can be readily combined with existing SPM toolboxes for estimating diffusion MRI parameter maps. From a user's perspective, the hMRI-toolbox is an efficient, robust and simple framework for investigating qMRI data in neuroscience and clinical research.


Introduction
Quantitative 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 (Cercignani et al., 2018;Assaf and Basser, 2005;Draganski et al., 2011;Lorio et al., 2014Lorio et al., , 2016aStüber et al., 2014;Callaghan et al., 2015a). 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 (Mugler and Brookeman, 1990) versus modified driven equilibrium Fourier transform, MDEFT (Deichmann et al., 2004) for T1-weighted anatomical images), sequence parameters (e.g. repetition time, TR, echo time, TE, 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, T 1 and T 2 , or the proton density, PD . qMRI accounts for these varied effects in order to increase the specificity of the estimated metrics and eventually quantify specific physical tissue properties (Cercignani et al., 2018;Lutti et al., 2010;Weiskopf et al., 2013). In qMRI, the estimated physical value has a direct meaning and is quantified in standardised units (e.g. T 1 in seconds) (Koenig et al., 1993). This standardised nature further increases the comparability across sites and time points (Deoni et al., 2008;Weiskopf et al., 2013), 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 (Callaghan et al., 2015a;Stikov et al., 2015) or a combination of qMRI with biophysical modelling (e.g. Henkelman et al. (2001), Assaf and Basser (2005), or Mohammadi et al. (2015)) 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, Weiskopf et al. (2015)).
The estimation of quantitative and semi-quantitative metrics commonly includes one or more of the effective transverse relaxation rate (R ⋆ 2 ¼ 1=T ⋆ 2 ), the longitudinal relaxation rate (R 1 ¼ 1=T 1 ), the proton density (PD), the magnetisation transfer (MT) saturation and a number of diffusion MRI (dMRI) metrics (Draganski et al., 2011). 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 (Sotiropoulos et al., 2013)) and processing software (see a summary in Soares et al. (2013)) are readily available for dMRI but less so for other qMRI techniques.
Consequently, 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 R ⋆ 2 , R 1 , PD and MT 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, R ⋆ 2 is often estimated using gradient recalled echo data acquired with multiple echo times (TE) to create high resolution maps that show strong contrast between different types of brain tissue (Bernstein et al., 2004). Similarly, PD and R 1 can be derived from two acquisitions varying the excitation flip angle (Wang et al., 1987;Deoni et al., 2005;Deoni, 2007;Schabel and Morrell, 2008;Helms et al., 2008aHelms et al., , 2011Liberman et al., 2013;Heule et al., 2015;Baudrexel et al., 2016). For accurate estimation of the latter qMRI parameters, an implementation must adequately correct for instrumental biases such as inhomogeneous transmit (Lutti et al., 2010) and receive fields (Volz et al., 2012;Mezer et al., 2016;Lorio et al., 2018). 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.
While 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 (Smith et al., 2006) and TRACULA (Yendiki et al., 2011) 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 (Dale and Sereno, 1993;Fischl et al., 1999;Fischl and Dale, 2000) to compare quantitative relaxation and susceptibility data on the cortex (Marques et al., 2017), or usage of the VBM framework (Ashburner and Friston, 2000) to process qMRI data across the whole brain on a voxel-by-voxel basis (Büchel et al. (2004) and Mohammadi et al. (2012) for dMRI metrics, and Table 1 for R ⋆ 2 , R 1 , PD and MT 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 Draganski et al. (2011) and developed for the comprehensive multi-parameter mapping (MPM) approach (Helms et al., 2008bWeiskopf et al., 2011Weiskopf et al., , 2013 to correct for potential error introduced by spatial smoothing. A 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 R ⋆ 2 , R 1 , PD, as well as semi-quantitative MT saturation maps. This framework enables timeefficient whole brain mapping with high isotropic resolution of 800 μm in 27 min (Callaghan et al., 2015b) or reduced MPM protocol (no MT saturation) at 1 mm isotropic resolution in 14 min at 3T (Papp et al., 2016) and has even enabled the acquisition of ultra-high-resolution quantitative maps with 400 μm resolution at 7T (Trampel et al., 2017). 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.
In this paper, we present the hMRI-toolbox, a comprehensive opensource toolbox that streamlines all the processing steps required to generate R ⋆ 2 , R 1 , PD and MT 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 R ⋆ 2 or variable flip angle data for mapping of R 1 and PD 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.

The MPM protocol
The MPM multi-echo protocol was introduced in Weiskopf and Helms (2008) and Weiskopf et al. (2013) for estimating the longitudinal relaxation rate R 1 , the effective transverse relaxation rate R ⋆ 2 , the proton density PD and the magnetisation transfer MT and generalises a number of acquisition protocols. It typically involves acquiring six to eight images at different echo times (TE) 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 R ⋆ 2 or variable flip angle data for mapping of R 1 and PD using multi-echo or single-echo data, comparable to e.g. DESPOT1 (Deoni et al., 2005). Example MPM acquisition protocols can be found at http://hmri.info.

Overview theory of MPM signal model
We 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.
The signal from the multi-echo PDw, T1w and MTw acquisitions can be described by the Ernst equation (Ernst and Anderson, 1966;Helms et al., 2008a, b). The effective transverse relaxation rate R ⋆ 2 can then be derived from the TE dependence of the signal. The unified description of the multi-echo data from all three contrasts into a single model, denoted Table 1 Review of studies related to the MPM acquisition protocol using predecessors of the hMRI-toolbox to make inference on myelin (My), iron (Fe), or the volume (Vol) measured by voxel-based morphometry (VBM). Note that all studies used Siemens MRI scanners. Abbreviations: R 1 ¼ longitudinal relaxation rate, R ð⋆Þ 2 ¼ (effective) transverse relaxation rate, PD ð*Þ ¼ (effective) proton density, MT ¼ magnetisation transfer saturation, C ¼ controls, P ¼ patients, dMRI¼ diffusion MRI, fMRI¼ functional MRI, MEG¼ magnetoencephalography, EEG¼ electroencephalography.

C
Review on cortical models for in vivo histology.
K. Tabelow et al. NeuroImage 194 (2019) 191-210 as ESTATICS , provides a more robust estimation of R ⋆ 2 with a higher signal-to-noise ratio compared to separate estimations (Fig. 1a). Using approximations of the signal equations for small repetition time TR and small flip angles α, the longitudinal relaxation rate R 1 , the apparent signal amplitude A Ã map (proportional to the proton density PD) and the magnetisation transfer MT can be estimated. At this point (Fig. 1b), the generated maps are biased by B 1 transmit f T (Fig. 1c) and receive f R (Fig. 1d) field inhomogeneities. The hMRI-toolbox provides correction methods for these bias fields based on specific B 1 transmit and receive field measurements or image processing methods. While f T influences the local flip angle and hence all three (R 1 , PD, MT) maps are affected, the RF sensitivity bias field f R only influences the PD map (in the absence of subject motion).
The toolbox can also handle the situation where only a subset of data is available. For example, R ⋆ 2 , R 1 and PD can still be estimated when no MTw acquisitions are acquired, R ⋆ 2 alone can be estimated when neither MTw nor T1w acquisitions are available (single multi-echo PDw data). R 1 , PD and MT 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 R 1 maps from other variable flip angle approaches, such as DESPOT1 (Deoni et al., 2005) or R ⋆ 2 maps from multi-echo data, such as certain susceptibility mapping/weighted imaging approaches.

Toolbox documentation and installation
The latest version of the toolbox can be downloaded from the hMRItoolbox page (http://hmri.info) as a zip file (containing the last official release) or by cloning the git repository (https://github.com/hMRIgroup/hMRI-toolbox) to keep up-to-date with the latest incremental developments.
Updated 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.
The 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.
The 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/.
In 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.

MPM example dataset
An 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 Callaghan et al. (submitted).

Organisation of the toolbox
The 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, Autoreorient 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.

Configure toolbox module
The 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).

DICOM import module
DICOM 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 hMRItoolbox implementation. Note that Philips-specific rescaling factor (Chenevert et al., 2014) 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, Fig. 1. Overview of qMRI map generation from the weighted imaging and reference MPM data. The signal S is modelled by the Ernst equation with an exponential decay depending on the echo time TE. The longitudinal relaxation rate R 1 , the effective transverse relaxation rate R ⋆ 2 , the proton density PD and the magnetisation transfer (MT) saturation are estimated from the data, using approximations for small repetition time TR and small flip angles α. The transmit and receive bias fields f T and f R are used to correct for instrumental biases. the metadata handling tools and BIDS compliance aspects (Gorgolewski et al., 2016) is provided in Appendix B.

Auto-reorient module
The 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 (Mazziotta et al., 1995(Mazziotta et al., , 2001a, is a common step that increases the consistency in individual head positions prior to normalisation or segmentation. For example, SPM's segmentation (Ashburner and Friston, 2005) 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.

Create hMRI maps module
The Create hMRI maps module computes quantitative as well as semiquantitative estimates of R ⋆ 2 , R 1 , PD and MT 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.

Multi-parameter input images
The 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 R ⋆ 2 , R 1 , PD and MT saturation (Helms et al., 2008a, b;Weiskopf et al., 2013Weiskopf et al., , 2014 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 R ⋆ 2 map. If T1w images are additionally provided, R 1 and PD maps will also be generated. The MT map estimation requires the acquisition of additional MTw images. The map creation also involves optional correction for B þ 1 (B 1 transmit) bias field f T (Section 3.7.2) and the B À 1 (RF receive sensitivity) bias field f R (Section 3.7.3) as well as spoiling imperfections (Yarnykh, 2010;Preibisch and Deichmann, 2009).

B 1 (transmit) bias correction
The map creation module includes the determination of B 1 transmit bias field maps (f T 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 (Lutti et al., 2010Weiskopf et al., 2013Weiskopf et al., , 2011Yarnykh, 2007;Chung et al., 2010).

Receiver RF sensitivity bias correction
Three options are available to correct for RF receive sensitivity bias (f R ) 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 intercontrast variation in RF sensitivity profile due to larger subject motion (Papp et al., 2016). Details on the different RF sensitivity bias correction methods can be found in Appendix C.4 and in the respective publications (Papp et al., 2016;Weiskopf et al., 2013).

Output
By default, the estimated quantitative maps are output into a Results subdirectory within the folder of the first PDw echo. Alternatively, a userdefined 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 Fig. 2. Left: After installation (Section 3.1), the hMRI toolbox can be started from the SPM menu of the batch editor. Five options include toolbox configuration, DICOM import, re-orientation of the data to MNI space, creation and processing of the hMRI maps. Right: Input mask for the map creation part of the toolbox. 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.

Process hMRI maps module
The 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 (Ashburner and Friston, 2005;Draganski et al., 2011;Lorio et al., 2016a), (2) non-linear spatial registration into common space (Ashburner (2007)) and (3) tissue-weighted smoothing (Draganski et al. (2011)), 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.

Statistical analysis
The 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., Draganski et al. (2011) andFreund et al. (2013). Additionally, the multiple parameter maps lend themselves to multi-variate analyses approaches as well (Draganski et al., 2011).

Discussion and outlook
This paper introduced the hMRI-toolbox, which is embedded in the  (Callaghan et al., submitted). The MPM protocol includes three multi-echo spoiled gradient echo scans with predominant T1-, PD-and MT-weighting achieved by an appropriate choice of the repetition time, the flip angle and the off-resonance MT pulse. Optional RF transmit (B þ 1 ) and receive (B À 1 ) field measurements can be added to the protocol, improving the quality of instrumental bias correction in the MPM maps. Alternatively, these reference measurements can be, to a limited extent, replaced by dedicated image processing steps that are provided by the toolbox. The map creation module generates PD, R 1 , MT and R ⋆ 2 maps. For each map, a JSON metadata file is created, which contains information about the processing pipeline of each image (see example in Table 3).

Table 2
Output files from the Create hMRI maps module using the SE/STE B 1 mapping and per-contrast RF sensitivity bias correction.

Results/Supplementary directory Description
hMRI_map_creation_rfsens_params.json RF sensitivity bias correction parameters (measured sensitivity maps) hMRI_map_creation_b1map_params.json B1 transmit map estimation: acquisition and processing parameters hMRI_map_creation_job_create_maps.json Create hMRI maps: acquisition and processing parameters hMRI_map_creation_mpm_params.json Acquisition and processing parameters used for the current job hMRI_map_creation_quality_assessment.json Quality assessment results (Appendix D) <firstSESTEfileName>_B1map. [niijjson] Estimated B1 bias field fT map (p.u.)

<firstSESTEfileName>_B1ref.[niijjson]
Anatomical reference for B1 bias field correction <firstPDfileName>_MTw_OLSfit_TEzero. [niijjson] MTw echoes extrapolated to TE ¼ 0 <firstPDfileName>_PDw_OLSfit_TEzero. [niijjson] PDw echoes extrapolated to TE ¼ 0 <firstPDfileName>_R2s. [niijjson] Estimated R ⋆ 2 map from simple exponential fit (PDw echoes) <firstPDfileName>_T1w_OLSfit_TEzero. [niijjson] T1w echoes extrapolated to TE ¼ 0 SPM framework and allows for the estimation and processing of four quantitative MRI parameter maps: the longitudinal and effective transverse relaxation rates R 1 and R ⋆ 2 , the proton density PD and the (semiquantitative) magnetisation transfer saturation MT. 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.
The 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).
Considerations 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 PD, R ⋆ 2 , R 1 and MT, 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.
At a given field strength, the R 1 contrast is determined by the microstructural 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 R 1 depends within limits linearly on these tissue properties (Fatouros and Marmarou, 1999;Fatouros et al., 1991;Kaneoke et al., 1987;Shuter et al., 1998;Donahue et al., 1994;Kamman et al., 1988;Fullerton et al., 1982;Gelman et al., 2001). The R ⋆ 2 and MT metrics can be used to describe the dependence of R 1 on these components (Callaghan et al., 2015a).
MT saturation is a proxy measure of the bound-pool water fraction (Helms et al., 2008b). 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. Freund et al. (2013) and Callaghan et al. (2014). Under equivalent conditions for the off-resonance pre-pulse the same MT 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 MT saturation measure as being semi-quantitative. The MT saturation map differs from the commonly used MT ratio (MTR; percent reduction in steady state signal) by explicitly removing the bias introduced by the spatially varying T 1 relaxation time and B 1 -transmit field (Helms et al., 2008b). Additional minor corrections for B 1 transmit field inhomogeneity in the MT maps were applied as described in Weiskopf et al. (2013). 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 (Callaghan et al., 2016b). Note that the MT 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. Battiston and Cercignani (2018)). A more direct measure of the bound-pool fraction is provided by quantitative MT (qMT), which requires more time-consuming data acquisition, typically limiting the qMT map's spatial resolution (e.g., 2 mm isotropic in Stikov et al. (2011)).
Sources 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, T ⋆ 2 -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.
The ideal imaging protocol includes dedicated measurements of transmit and receive field inhomogeneities, which are typically based on Fig. 4. Overview of the spatial processing module. It consists of three steps: (1) segmentation, (2) highly parametrised nonlinear spatial registration, and (3) tissueweighted smoothing. The segmentation step (1) uses novel tissue-probability maps (TPMs), designed to take advantage of the better contrast of the MPMs for improved segmentation (exemplified for the deep-grey matter by the arrow in (1)). The non-linear spatial registration step into common space (2) reduces inter-individual anatomical differences (exemplified for R 1 maps of subject one and two, S1 and S2, respectively). To further reduce residual anatomical differences (see magnification boxes in (2)) and enhance statistical inference, the qMRI maps can be spatially smoothed (3) using the voxel-based quantification (VBQ) smoothing procedure. As compared to Gaussian smoothing, VBQ smoothing avoids bias in the qMRI maps (e.g. see arrows in (3), highlighting a rapid decline in R 1 values at tissue boundaries only after Gaussian smoothing. The VBQ smoothing is detailed in Eq. (C.3) and Draganski et al. (2011)). The sub-figure (1) has been adapted from (Lorio et al., 2016a), sub-figures (2) and (3) from Mohammadi and Callaghan (2018).
customised sequences (Fig. 5a-d for RF receive bias correction, and Fig. 6a for B 1 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 (Lutti et al., 2010 or 3D_AFI (Yarnykh, 2007) for correction of B 1 transmit effects can be replaced by manufacturer's service sequences such as TFL_B1_map (Chung et al., 2010) or RF_map (as examples on Siemens scanners). The toolbox provides the option to process data from several different B 1 transmit bias field mapping techniques or to use B 1 transmit bias field maps pre-computed outside the toolbox (see section 3.7.2 and Appendix C.3).
The 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 betweencontrast 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 (Hoult, 2000). 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, B 1 transmit field mapping techniques can be inaccurate. For a comparison of the frequently used B 1 -transmit field mapping techniques and a description of their respective accuracy and sources of uncertainty we refer to Lutti et al. (2010) and Pohmann and Scheffler (2013).
When 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 (Ashburner and Friston, 2005) for B 1 transmit bias correction (UNICORT, Weiskopf et al., 2011;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 TR (Helms et al., 2008a). This requires no additional acquisition time but produces quantitative maps of lesser accuracy with some residual receive and/or transmit field modulation Baudrexel et al., 2016) compared to a correction with measured references. Note that the Unified Segmentation approach, whether applied for B 1 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 Weiskopf et al. (2011) for UNICORT optimisation). For more details, see (Callaghan et al., 2016a) and the hMRI wiki (the latter also provides further information on customized usage).
The proposed MPM protocol uses RF and gradient spoiling to minimise undesired transverse net magnetisation. Imperfect spoiling, which PDw and T1w images corrected by a single sensitivity map, acquired directly before the PDw images (b: ΔPD S1 ), the MTw images (c: ΔPD S2 ) and the T1w images (d: ΔPD S3 ), respectively. Due to the large overt movement preceding the acquisition of the MTw images and corresponding RF sensitivity measurement (see Callaghan et al. (submitted)), there is a large discrepancy between head position for that specific RF sensitivity measurement on the one hand and head positions for the PDw and T1w images used to generate the PD map on the other hand. As a result, errors in (c) are much larger than in depends on the precise sequence protocol settings, can leave a residual bias in the R 1 map if no further correction is used (Preibisch and Deichmann, 2009;Yarnykh, 2010). 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).
The estimation of PD maps can be biased by T ⋆ 2 relaxation effects if not accounted for (Eq. (A.4) in Appendix A). Two correction methods, based on extrapolation of the data to TE ¼ 0 (Ellerbrock and Mohammadi, 2018) or relying on the estimated R ⋆ 2 maps respectively, are implemented in the toolbox (Balteau et al., 2018). These correction methods require a number of echoes to be acquired for a robust fit of the exponential decay to derive TE ¼ 0 magnitude images and R ⋆ 2 maps. In the case of a single-echo variable flip angle dataset, the T ⋆ 2 -weighting correction cannot be applied and the estimated PD and (to a much lesser extent) the R 1 and MT maps will be biased by T ⋆ 2 -modulations. Also, R ⋆ 2 maps are biased in areas with severe susceptibility artifacts (Yablonskiy, 1998). Note that the JSON metadata file associated with the respective PD parameter map contains information about the processing steps and thus of potential T ⋆ 2 modulation. Head 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 (Callaghan et al., 2015b;Weiskopf et al., 2014;Reuter et al., 2015). 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 (Rosen et al., 2018). The hMRI-toolbox provides summary measures of head motion within and between the acquisitions of each image volume (intra-and inter-scan motion) (Castella et al. (2018) 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 (Castella et al., 2018). 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.
Spatial 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 (Ashburner and Friston, 2000;Ridgway et al., 2008;Focke et al., 2011). 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.
Residual 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 (Ashburner, 2007). 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.
The choice of the appropriate smoothing kernel and its performance compared to alternative methods (e.g. TBSS (Smith et al., 2006) or TSPOON (Lee et al., 2009)) is still a subject of active research. An alternative method for reducing spatial misregistration in the cortex might be surface-based registration algorithms (Davatzikos and Bryan, 1996;Drury et al., 1996;Thompson and Toga, 1996;Fischl et al., 1999). However, a comparison study by Klein et al. (2010) 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 (Canna et al., 2018).
hMRI-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 (Lee et al., 2017(Lee et al., , 2018, and the MPM framework together with the hMRI-toolbox will be used in a multi-site clinical trial (NISCI trial, (Seif et al., 2018),) including Philips and Siemens MR systems.
Multi-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 MT saturation map, which will require appropriate interscanner calibration (Volz et al., 2010;Seif et al., 2018). Moreover, not every DICOM to NIfTI conversion software appropriately handles image intensity scaling, as reported e.g. for Philips data (Chenevert et al., 2014), leading to spurious intensity differences affecting the quantification. Making the hMRI-toolbox data formats fully BIDS compliant (Gorgolewski et al., 2016) by defining the ontology of acquisition parameters  Fig. 5f. 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.
Applicability 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. laminarspecific information (see 400 μm new generation MPM, Trampel et al. (2017)), while also posing new challenges (e.g. B þ 1 field inhomogeneities) that the hMRI-toolbox will have to address.
Future 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 hMRItoolbox for different MRI scanner platforms above).
Sensitivity 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 Callaghan et al., 2015a;Castella et al., 2018). Retrospective robust estimation of R ⋆ 2 parameters (outlier rejection) has been suggested in Weiskopf et al. (2014) 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 Tabelow et al., 2017) are important future improvements of the hMRI-toolbox.
Currently, 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.
Extensions 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 (Acosta-Cabronero and Callaghan, 2017;Acosta-Cabronero et al., 2018), 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 (Stikov et al., 2015;Mohammadi et al., 2015;Ellerbrock and Mohammadi, 2018). 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. MT or PD maps) with the axonal-sensitive maps obtained with existing SPM tools (e.g. the ACID toolbox (Mohammadi et al., 2012;Tabelow et al., 2015;Ruthotto et al., 2013;Mohammadi et al., 2014) and the DTI & Fiber tools (Reisert et al., 2013)).

Conclusion
The 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. For a single excitation and echo time TE, the signal intensity S is given by the Ernst equation for an ideal spoiled gradient echo acquisition (Ernst and Anderson, 1966;Helms et al., 2008a):

Acknowledgements
We denote the repetition time by TR and the flip angle by α. A is proportional to the equilibrium magnetisation and thus to the proton density PD (Helms et al., 2008a).
For fast estimation of the R 1 and A parameters from (A.1) and a single echo of a T1w and a PDw acquisition, or the mean across multiple echoes for increased SNR, Helms et al. (2008a) introduced a linearisation of its exponential terms for short TR (i.e. R 1 Á TR ≪ 1): leading to, after averaging over echos, When S is the averaged signal over several echoes, the sum in Eq. (A.4) is calculated across all echoes included in the averaged signal. This introduces an R ⋆ 2 bias which increases with the number of echoes and depends on the TE values included in the average. The R ⋆ 2 bias in A ⋆ can be evaluated using the R ⋆ 2 estimate and corrected for. Alternatively, S can be the extrapolated signal at TE ¼ 0, which is free from R ⋆ 2 weighting. Both the echo averaging with R ⋆ 2 correction method and the TE ¼ 0 extrapolation method are implemented and available in the toolbox, and perform similarly (Balteau et al., 2018). The toolbox extrapolates to a TE of 0 ms by default.
Using the approximation sinðαÞ $ α and cosðαÞ $ 1 À α 2 =2 for small flip angles and neglecting the term α 2 Á R 1 Á TR, Eq. (A.1) simplifies to: from which R 1 and A ⋆ are estimated using the signals from the PDw and T1w sequences: where signals S T1 and S PD are either averaged over a number of echoes (Helms et al., 2008a), or extrapolated to TE ¼ 0, as explained above. In the first case, a P TE e ÀR ⋆ 2 Á TE correction factor must be calculated to derive A (Balteau et al., 2018). A ¼ A ⋆ otherwise meaning that the R ⋆ 2 bias remains.

Appendix A.2. Estimation of the magnetisation transfer
The MT (magnetisation transfer) saturation measure is derived in analogy to a dual-excitation spoiled gradient echo sequence where the second excitation is replaced by the MT pulse. Using a heuristic approximation for small read-out flip angles α MT , the closed form solution for the acquired signal S can be written (Helms and Hagberg, 2009;Helms et al., 2008b) as where TR MT ¼ TR 1 þ TR 2 is the total repetition time, and δ describes the MT saturation (more generally denoted MT in this paper). From an additional MTw signal S MT we obtain, using again the linearisation of the exponential term for short TR MT (i.e. R 1 Á TR MT ≪ 1) and the small angle approximation for the trigonometric functions: see Helms et al. (2008b). Note that the parameter maps R 1 , R ⋆ 2 , A and MT can alternatively be estimated from the ESTATICS model using analytic formulas and not relying on the linearisation outlined above . with indicator variables I PD , I T1 , and I MT for the differently weighted acquisitions. The signal intensities, e.g. S PD , at TE ¼ 0 are given by with the flip angle and repetition time of the PDw sequence and correspondingly for S T1 and S MT . From Eq. (A.10) R ⋆ 2 can be estimated using all echoes. Appendix A.4. Corrections for transmit and receive field inhomogeneities Due to inhomogeneities of the B 1 transmit field the local flip angle deviates from its nominal value (Helms et al., 2008a). Using a correction field f T , the effective local flip angle is given by For details on the correction field f T estimation see Section 3.7.2 and C.3.
For the MTw acquisition we denote by α sat the off-resonance Gaussian shaped RF pulse used for MT weighting. The B 1 transmit bias field f T ðxÞ leads to a local correction of α sat , too, and thus to a residual bias field for the estimated δ even though δ is largely independent of the local transmit field bias unlike other MT metrics such as MTR (Helms et al., 2008b;Callaghan et al., 2014). In Helms (2015) the heuristic correction factor for δ was found to be with B Á α sat ¼ 0:4 for the typically used value α sat ¼ 220 in MPM sequences .
Finally, the multiplicative RF sensitivity bias field f R (see Section 3.7.3 and Appendix C.4) for the correction of the signal amplitude maps A: has to be measured or estimated from the data (for details see Section 3.7.3). The A ⋆ maps can then be normalised to proton density PD maps by calibration. To this end, the mean PD value in white matter is assumed to be 69 percent units (p.u.) (Filo and Mezer, 2018). Appendix A.5. Spoiling corrections Further, imperfect spoiling requires an additional correction when estimating R 1 (Preibisch and Deichmann, 2009), to account for deviation from Ernst equation (A.1). The correction is performed using the following expression where P a ðf T Þ and P b ðf T Þ are quadratic functions in f T . Their coefficients depend on the specific sequence and can be determined by simulation (Preibisch and Deichmann, 2009). Further refinements of this correction include diffusion effects on signal spoiling (Yarnykh, 2007;Lutti and Weiskopf, 2013;Callaghan et al., 2015c).

Appendix B. DICOM import and JSON metadata
To support the flexibility of the hMRI processing pipeline and the traceability of the data at every step, a number of parameters from acquisition to statistical results must be retrieved and stored. Therefore, DICOM import and JSON metadata handling functionalities have been implemented with the following objectives: 1. Retrieving the full DICOM header as a readable and searchable Matlab structure (metadata), 2. During DICOM to NIfTI conversion, storing the metadata structure as separate JSON file (alongside the NIfTI image file) and/or as extended header into the NIfTI image file, 3. Handling metadata in order to set, get and search parameters.
With the current implementation, these objectives are achieved: the hMRI-toolbox handles and stores data acquisition and processing parameters as JSON-encoded metadata alongside brain imaging data sets. The implementation rests on the SPM12 implementation of DICOM tools (spm_di-com_header(s), spm_dicom_essentials, spm_dicom_convert) and relies on spm_jsonread and spm_jsonwrite for handling JSONencoded metadata.
From SPM12 version r7219 (released 16 November 2017) and in later versions, functionalities described in points 1. and 2. above (excluding the extended header option though) have been integrated to the DICOM import module in SPM12 (implementation from the hMRI-toolbox). Note that some features of the hMRI-toolbox metadata library (for metadata handling, in particular to set, get and search the metadata) are not available in the SPM12 distribution.
Detailed syntax and many examples are provided in the hMRI-toolbox Wiki, we briefly describe the main features and the metadata structure below. Appendix B.1. DICOM to NIfTI conversion.
The spm_dicom_convert script has been modified to achieve the above goals. An additional input argument has been set to deal with the new JSON metadata options. If JSON metadata storage is enabled, DICOM header information is stored as a JSON-formatted structure either as a separate JSON file or as extended NIfTI header (or both). If omitted or disabled, the DICOM to NIfTI conversion proceeds using exactly the same implementation as in the previously mentioned SPM12 distribution. The output NIfTI images, extended or not, have a valid NIfTI-1 format compatible with standard NIfTI readers (MRICron, SPM, FreeSurfer). Note that FSL up to version v5.0.9 does not support extended headers. For that reason, the extended header option, although suggested in (Gorgolewski et al., 2016) is currently not recommended and disabled by default. Instead, a separate JSON file is created. Appendix B.2. DICOM to NIfTI conversion -SPM Batch.
The DICOM to NIfTI conversion with the above features is available through the hMRI-toolbox in the SPM Batch GUI (SPM > Tools > hMRI Tools > DICOM Import). The new JSON metadata format field provides the user with the following options: None Separate JSON file (default) Extended nii header Both Appendix B.3. The metadata structure.
Metadata are simple MATLAB structures, hence quite flexible and modular. The structure of the metadata is divided into the following two main fields: Acquisition parameters (hdr.acqpar): this branch contains the original DICOM header. It should be kept unchanged and is created when importing DICOM images to NIfTI. It is dropped for processed images relying on several input images. History (hdr.history): is a nested structure containing procstep: a structure describing the current processing step and containing: * descrip: a brief description of the processing applied (e.g. DICOM to NIfTI conversion, realign, create map) * version: version number of the processing applied, for traceability * procpar: processing parameters (relevant parameters used to process the datae.g. for SPM processing, all the settings included in the MATLAB batch job except for the input images. input. An array listing the input images used for the processing. Each input (hdr.history.input(i)) is a structure containing: * filename: the filename of the input image * history: the history structure from the input image output: a structure containing * imtype: image type of the current output (e.g. R 1 , FA, MT, ADC) including a summary description of the processing steps (see TableB 3 The JSON transcription of the MATLAB structure is easily readable and searchable by opening the extended NIfTI or separate JSON file with a text editor. Furthermore, metadata set, get and search tools have been implemented. These tools have been added into the spm12/metadata directory. Typing help <function> in MATLAB gives detailed syntax. Below is a list of the most useful scripts and their summary descriptions. For a complete description of the implemented tools, please refer to the hMRI-toolbox Wiki pages. get_metadata: returns the MATLAB structure described above. set_metadata: to write (insert, modify, overwrite) metadata in JSON files and extended NIfTI headers. In general, metadata are initialised during the DICOM to NIfTI conversion and further modified as the processing progresses.
get_metadata_val: returns pairs of structure fields and values agreeing with search criteria. find_field_name: recursively searches a MATLAB structure for a specific field name (can be applied to any MATLAB structure). Appendix B.5. General usage of the metadata and BIDS compatibility The DICOM import with metadata storage is available in the SPM distribution since release r7219 (November 2017) allowing for a wide use of such metadata. It provides most of the parameters required for map creation and processing. The hMRI-toolbox also provides metadata handling functionalities to retrieve parameter values and store processing parameters in the JSON-encoded metadata file. The format used for the metadata provides flexibility for using a wide range of protocols, and traceability of every parameter used to generate the results (see for example Table 3). Table 3 JSON Metadata example: Metadata are stored as a Matlab structure in the JSON file associated to each output (NIfTI) image. The metadata structure can easily be retrieved using the following SPM/Matlab command: metadata ¼ spm_jsonread(<name of the JSON file>). The structure includes processing parameters, dependencies of the result on other input images as well as a summary description of each output image. As an example, the image type of each output R 1 image used to generate Fig. 6 is provided in the Table. The highlighted lines (bold font) describe the B 1 transmit bias correction, illustrating the user-friendly traceability of the processing steps used to generate an image. The same information can be read directly from the JSON file using a text editor. Outside the hMRI-toolbox, metadata and metadata handling tools implemented in the hMRI-toolbox can prove useful for many purposes. For quality control, metadata can help checking for protocol consistency by comparing essential parameters across sessions, subjects and sites. When processing data acquired with several protocols on several scanners from different vendors (multi-centre studies), metadata provide the information to automatically extract acquisition parameters to be used as regressors of no interest in a statistical analysis (e.g. MR scanner, field strength, head coil, …). Metadata also facilitate the sorting of the data between different data types (e.g. functional MRI, structural MRI, diffusion images, etc.) together with retrieving BIDS essential parameters (e.g. repetition time, echo time, flip angle, session and series numbers, etc.), therefore providing all the information required to make the data archiving and storage fully BIDS-compatible.
A BIDS extension proposal (BEP) is currently under discussion, that will define the parameters necessary to describe structural acquisitions with multiple contrasts (see https://bids.neuroimaging.io/for updates on BIDS specifications and extension proposals), and bring the metadata structure of such data in line with the BIDS recommendations (Gorgolewski et al., 2016). The extension will include all acquisition parameters required to fully describe and process data acquired with the MPM protocols, in a standard, platform-and vendor-independent way. With such BIDS extension, the hMRI-toolbox will be adapted and made fully BIDS-compliant, from the structure of the JSON-encoded metadata to the structure of the output data. The required adaptation will be achieved smoothly thanks to the metadata handling functionalities, without affecting the processing part of the hMRI-toolbox, and enabling the implementation of fully automated tools retrieving the BIDS-organised data, generating the maps and applying spatial processing.
Appendix C. Detailed description of the hMRI modules Appendix C.1. Configure toolbox module The Configure toolbox module must be run before any other module for the default parameters to take effect. The default settings are defined and described in the config directory. The default files therein should never be modified. Template files for customised default parameters are stored in the config/local directory and can be modified and saved with meaningful names by the user. Note that most parameters should only be modified by advanced users who are aware of the underlying theory and implementation of the toolbox.
For brevity, a full description of each customisable setting is not provided here. A detailed description of each parameter may be found in the hmri_local_defaults.m help and examples are provided in the toolbox Wiki pages. Note that the default settings for B 1 transmit bias correction (also including acquisition and processing parameters) can be customised separately within the Create hMRI maps module (Appendix C.3). Their customisation must follow the same guidelines as provided in this section.
While settings need to be defined either in the standard or customised default files, acquisition parameters can be retrieved from the input images if included in the metadata. When no metadata are available or not yet properly handled by the toolbox (e.g. customised sequences), acquisition parameters specified in the default files are used as a fallback solution. The use of metadata is strongly recommended (see Appendix B), to automatically retrieve acquisition parameters and to avoid incorrect processing.

Appendix C.2. Auto-reorient module
This module performs a rigid-body alignment of a subject's structural image into the MNI space, applying the estimated pose to a set of other images. The code makes use of the spm_affreg function and canonical images available in SPM12. The user must provide one structural image from the subject and a corresponding canonical image for co-registration, and a series of other images that need to be kept in alignment with the structural image (typically, all the images acquired during a given imaging session). The structural image should have good signal-to-noise ratio with a high white matter/grey matter contrast-to-noise to ensure a robust and reliable co-registration. The SPM canonical images can be selected from the SPM/canonical directory or any other user-defined template (e.g. atypical population-driven template), provided it is already oriented according to MNI. Since the registration is based on matching similar intensities, it is important that the contrast of the reference image be close to that of the canonical image. The following list is provided as guidelines: If no specific output directory is selected, the original orientation of the images will be overwritten. If an output directory is selected then the images are copied there before being reoriented, therefore leaving the original data untouched. The reoriented images can be collected as dependencies for further processing in the batch interface, either all at once (Dependencies > Grouped) or as individual output images (Dependencies > Individual), which is more convenient to connect the Auto-reorient module to the next Create hMRI maps module. Appendix C.3. B 1 transmit field corrections Within the Create hMRI maps module, several methods for B 1 transmit field correction are implemented. Depending on the choice of the specific method the GUI requires the user to provide adequate input files, which are described below. In order to proceed to the B 1 map calculation, number of acquisition parameters must also be either retrieved from the metadata or provided to the toolbox via a standard or customised default B 1 file (see below and in Appendix C.1).
The sequence uses a set of pulses with nominal flip angles α, 2 Á α and α to create pairs of spin echo and stimulated echo images. All consecutive pairs of SE/STE images corresponding to the different nominal flip angle values α of the SE/STE RF pulse must be loaded as B1 input. B 0 field mapping images must also be provided for correcting distortions in the EPI images. Both magnitude images and the pre-subtracted phase image must be selected, in that order, as B0 input.
A pair of magnitude images acquired with two different repetition times (TR) must be loaded as B1 input. The hMRI-toolbox then calculates the B 1 transmit bias field map as described in Yarnykh (2007).

tfl_b1_map
TFL B 1 mapping (Chung et al., 2010). For this method, the batch interface requests a pair of images (one anatomical image and one flip angle map, in that order) from a service sequence by Siemens (version available from VE11 on) based on a turbo flash (TFL) sequence with and without a pre-saturation pulse (Chung et al., 2010). The flip angle map used as input contains the measured flip angle multiplied by 10. After rescaling (p.u.) and smoothing, the output f T map is ready to be used for B 1 transmit bias correction.

Rf_map
This requires a pair of images (one anatomical image and one pre-processed B 1 map, in that order) from a service sequence by Siemens (Erlangen, Germany). These are based on the acquisition of a spin-echo/stimulated echo as used for method 3D_EPI above (Lutti et al., 2009. Rescaling and smoothing is applied to the input B 1 map (a.u.) to generate a f T map (p.u.) suitable for B 1 transmit bias correction.

Pre-processed B1
Any B 1 transmit bias field map pre-calculated using one of the above methods or another method can be used as pre-processed B1 input. The user must select one anatomical reference (for co-registration) and one B 1 map (in p.u.), in that order. The anatomical reference is typically a magnitude image from the B 1 mapping dataset, and must be in the same space as the B 1 transmit bias field map.

UNICORT
Data driven estimation of the B 1 transmit bias field map using the UNICORT approach . If none of the above mentioned methods applies because no appropriate B 1 transmit bias field mapping data were acquired, the UNICORT option is recommended to remove the transmit field bias in the R 1 map and estimate the B 1 transmit bias field map from the qMRI data. Optimal UNICORT processing parameters are RF transmit coil dependant. The default parameters have been optimised for Siemens 3T TIM Trio scanners and may need specific optimisation for other RF transmit coils. See the hMRI-toolbox Wiki for guidelines.

No_B1_correction
With this option, no correction for transmit inhomogeneities is applied, i.e. f T 1. For each B 1 mapping protocol (1-4), one *_B1map (in p.u.) and one *_B1ref (for anatomical reference) file is created and saved in the Results/ Supplementary directory (see Table 2). The output images are associated with JSON metadata including the type of B 1 map processed and the setting used. For UNICORT, only the *_B1map in p.u. Is created (in the same space as the qMRI maps) with no anatomical reference. For pre-processed B 1 transmit bias field maps provided by the user, the input files are not renamed to match the naming convention above. Default acquisition and processing parameters for the UNICORT, 3D_EPI and 3D_AFI methods can be customised at this point by providing the Create hMRI maps module with a customised configuration file (see Appendix C.1 and the hMRI Wiki pages for guidelines). Except for UNICORT, where the B 1 transmit bias field is directly derived from the qMRI data, and the no_B1_correction option, where no anatomical reference is given, the hMRI-toolbox co-registers the anatomical reference image with the multi-echo spoiled gradient echo data before applying the B 1 transmit bias correction. Appendix C.4. Receive (RF) sensitivity bias correction In this paragraph, the RF sensitivity bias correction methods based on measured RF sensitivity maps (Single or Per contrast options) as well as the data driven method (Unified Segmentation option: no input sensitivity map required) will be explained in detail.
Receive field sensitivity measurements. Since the receive field sensitivity of the body coil is assumed to be flat over the spatial extent of the head (Pruessmann et al., 1999), if the same anatomy is imaged with the head coil and the body coil sequentially, using the same acquisition parameters and assuming no motion, then the ratio of these two scans is proportional to the net head coil receive sensitivity field f R : where S HC is the signal acquired with the head coil, f R is the receive sensitivity field of the head coil, S BC is the same signal acquired by the body coil, C BC is the receive sensitivity field of the body coil assumed to be approximately constant, and S 0 is the signal specific to the underlying anatomy and acquisition parameters. If receive field sensitivity measurements are available for each imaging contrast (one receive sensitivity field acquired either before or after each of the PDw, T1w and MTw contrasts), inter-scan variation of the sensitivity modulation can be accounted for and used to optimally correct for the combined inter-scan motion and sensitivity modulation effects in the quantitative maps (Papp et al. (2016) and Fig. 5a). Therefore, the implemented RF sensitivity correction method (Per contrast option) combines a correction for motion-related relative receive sensitivity variations with rigid-body realignment.
The low resolution measurements from the head and body coils (in this order) for each T1w, PDw and MTw acquisition are expected as inputs to the RF sensitivity bias correction. A receive field sensitivity map is calculated according to Eqs. C.1 and C.2 for each contrast, and applied to the corresponding contrast to correct for any coil sensitivity driven signal intensity modulation.
If RF sensitivity measurements are missing for some of the contrasts, or if only one single measurement is available, the same RF sensitivity bias correction can still be applied to all contrasts by choosing the corresponding (Single) option in the GUI. Only one pair of head coil and body coil images is required. Receiver field inhomogeneities will be corrected for although not accounting for inter-scan motion (Fig. 5b-d). Note that the measured RF sensitivity maps are modulated by the receive field of the body coil, which is not corrected for (Fig. 5e).
Data driven receive field estimation. If no receive field sensitivity map has been acquired, the Unified Segmentation option can be selected in the GUI. The receive field bias is then estimated from the calculated A map using SPM unified segmentation (US) approach (Ashburner and Friston, 2005). Accurate bias field estimation requires the use of a brain mask combining the white and grey matter probability maps. The latter are derived by segmentation of the calculated MT map (because of its higher contrast in the basal ganglia, see  or alternatively (if no MTw images are available), by segmentation of the calculated R 1 map. The optimal US parameters depend on the RF receive coil and its sensitivity profile. Fine adjustment of the US parameters might be required for a specific RF receive coil. Appendix C.5. Process hMRI maps module This section describes the three main operational steps of the spatial processing pipeline for hMRI data: segmentation, diffeomorphic registration and tissue-weighted smoothing. Furthermore, it provides details on the fully integrated processing pipeline, which is proposed in this toolbox as a separate sub-module to facilitate standard spatial data processing without the need to combine the individual steps.
The processing tools and user interface can be used to process parameter maps generated from a series of subjects, assuming that each subject has the same small number of images (for example the four MT saturation, PD, R 1 , and R ⋆ 2 maps). Therefore each module takes as input several series of images sorted per image type (for example the MT from all the subjects). The images in each series must follow the same subject-based ordering.
Appendix C.5.1. Segmentation This step relies on the latest unified segmentation (US) algorithm as available in SPM12 (Ashburner and Friston, 2005) and is interfaced through a single module. Compared to the SPM12-Segment module, the hMRI-Segment module allows the successive processing of a series of subjects.
The segmentation itself relies on one structural image per subject, typically an MT map, R 1 map or T1w image, and a set of tissue probability maps (TPMs). The TPMs used by default in the hMRI-toolbox were specifically derived from multi-parametric maps (Lorio et al., 2014) denoted as extended TPMs (eTPM).
By default (adjustable in the batch interface), the segmentation step generates a series of tissue class images in native space: the c1/c2/c3 images -for gray matter (GM), white matter (WM) and CSF, as well as Dartel-imported rc1/rc2 images. The segmentation module also outputs the GM and WM tissue class images warped into MNI space (classic elastic deformation, i.e. not with Dartel) with (wmc1/wmc2 images) and without (wc1/wc2 images) Jacobian modulation. Finally, the estimated segmentation and deformation fields are available as a.mat file and 4D y_*.nii files respectively. See the upper part of Table 4 for a complete file list and description of the output.

Table 4
Output files from the Process hMRI maps modules. The file names are based on the file name of the map or image used as input to the Segmentation, Dartel > Run Dartel (create Templates), Dartel > Normalise to MNI space, and Smoothing steps. Note that wc*<InputFileName> [niijjson]. images can be obtained either with the simple warp obtained with the Unified Segmentation or following Dartel with the Dartel > Normalise to MNI space module.

mwc1<segmInputFileName>.[niijjson]
Tissue class 1 warped into MNI space, with the simple warp obtained with the US, and modulated by the determinant of the Jacobian, i.e. accounting for local change of volume. Prefixes mwc1/mwc2 correspond to GM/WM respectively. y_<segmInputFileName>.nii Deformation field, i.e. warps, obtained from US <segmInputFileName>_seg8.mat Segmentation and warping parameters, obtained from US Diffeomorphic registration (Proc. hMRI -> Dartel -> Run Dartel (create Templates)) u_rc1*_Template.nii Flow field image, one per subject, estimated by Dartel from the rc1/rc2 images. Template_*.nii 7 template images, numbered from 0 to 6, created by Dartel from the rc1/rc2 images of all the subjects. Diffeomorphic registration (Proc. hMRI -> Dartel -> Normalise to MNI space) w*<normInputFileName>. [niijjson] Image warped into MNI space following Dartel, using the estimated flow field and an affine transformation. This would be typically a qMRI map that should not be modulated to account for volume changes. mw*<normInputFileName>. [niijjson] Image warped into MNI space and modulated by the determinant of the Jacobian, i.e. accounting for local change of volume. This would typically be any tissue probability map (i.e. image with a measure whose total amount over the brain volume should be preserved) to be used after smoothing for a VBM analysis. Tissue-weighted smoothing (Proc. hMRI -> Smoothing) wap1_<smooFileName>. [niijjson] Tissue-weighted smoothing for tissue class 1. s*<FileName>. [niijjson] (Any) image smoothed with a standard Gaussian filter (here for comparison, not in hMRI). alignment with MNI space. Unlike to the Dartel toolbox, no smoothing is applied at this stage, but a specific smoothing operation is applied at the next Tissue-weighted smoothing step. See the third part of Table 4 for the file list and description. The machinery behind those modules are exactly the same as the 'Dartel toolbox' available elsewhere in SPM. They have only been included in the hMRI-toolbox for ease of use. Moreover the user-interface of the Normalize to MNI module was slightly extended to more efficiently handle the multiple parameter maps and tissue class images of each subject. Notably, in contrast to many typical VBM applications, quantitative parameter values are not modulated (to account for volume changes) during this normalization.
These parameters are population-and protocol-dependent which is why they are referred to as relative indicators.