Fractional order magnetic resonance fingerprinting in the human cerebral cortex

Mathematical models are becoming increasingly important in magnetic resonance imaging (MRI), as they provide a mechanistic approach for making a link between tissue microstructure and signals acquired using the medical imaging instrument. The Bloch equations, which describes spin and relaxation in a magnetic field, is a set of integer order differential equations with a solution exhibiting mono-exponential behaviour in time. Parameters of the model may be estimated using a non-linear solver, or by creating a dictionary of model parameters from which MRI signals are simulated and then matched with experiment. We have previously shown the potential efficacy of a magnetic resonance fingerprinting (MRF) approach, i.e. dictionary matching based on the classical Bloch equations, for parcellating the human cerebral cortex. However, this classical model is unable to describe in full the mm-scale MRI signal generated based on an heterogenous and complex tissue micro-environment. The time-fractional order Bloch equations has been shown to provide, as a function of time, a good fit of brain MRI signals. We replaced the integer order Bloch equations with the previously reported time-fractional counterpart within the MRF framework and performed experiments to parcellate human gray matter, which is cortical brain tissue with different cyto-architecture at different spatial locations. Our findings suggest that the time-fractional order parameters, {\alpha} and {\beta}, potentially associate with the effect of interareal architectonic variability, hypothetically leading to more accurate cortical parcellation.


Introduction
Magnetic resonance imaging (MRI) is a routinely used medical imaging modality known for its exquisite soft-tissue contrast. The power of this imaging modality arises from its sensitivity to changes in tissue composition, interpreted as changes in texture, intensity, shape and size of structures within images. Knowledge of the drivers behind the changes allows accurate diagnosis, treatment planning and monitoring of patients presenting with a range of diseases and disorders. In the brain, MRI-based image contrast is primarily due to microstructural differences in gray and white matter tissue. As such, assessment of diseases and disorders relates to how additional information appears in images (e.g. tumours), or how gray-white matter shape and size, topological features and spatial intensity vary with disease.
The current trend in MRI has been to supplement the qualitative image repertoire with quantitative, biologically relevant maps of tissue properties to improve diagnostic, prognostic and treatment planning accuracy and specificity. This necessitates the research and development of analytical methods capable of inferring microstructural information from mm-scale measurements to bypass resolution limits imposed by the MRI hardware. The need for reproducible, quantitative MRI has prompted the development of methods capable of linking mm-scale MRI measurements with biologically interpretable information [1][2][3][4][5][6][7]. Here, an appropriate mathematical model in conjunction with specifically collected MRI data lead to model parameters sensitive at the microscale. MRI voxel factors including myelin and iron content, tissue density, composition and microstructure orientation have been estimated in this way, e.g. [1,3,[8][9][10].
Classically, the Bloch equations describes how the MRI signal evolves over time as a function of spin-lattice ( 1 ) and spin-spin ( 2 ) relaxation. The equation holds for the case of an isotropic material, for example a sample made up of water alone. However, when relaxation occurs in a complex tissue structure with tissue comprising of multiple constituents, the integer order differential equations describing sample magnetisation with a mono-exponential signal solution is no longer able to capture the trend in the signal as a function of time [11]. In this instance, the temporal MRI signal is said to deviate away from the expected mono-exponential trend. Models involving multiple exponentials (e.g. [12][13][14]) and time-fractional order derivative representations (e.g. [11,15]) have been proposed to better explain trends in MRI signals generated in real tissues. The reader is referred to a recent review article for a comprehensive overview of anomalous relaxation processes in MRI [16].
When using mathematical models applied to MRI data, the estimated model parameters are assumed to incorporate information on tissue microstructure. As such, a model needs to either exist or be developed, and then fitted to specifically collected data. Estimation of model parameters involves a non-linear solver, and often a good initial guess has to be made for the parameters to converge to realistic values [11]. This approach has shortcomings, including the potential for overfitting and making sure the non-linear solver converges to a realistic solution. A different approach is to generate a socalled dictionary based on discrete parameter choices. For each set of model parameters, a simulated signal based on the Bloch equations can be generated and matched with the observed signal. This latter framework is referred to as magnetic resonance fingerprinting (MRF) [17], and has become a routine tool in MRI. The parameters associated with the best matched signal are said to be those most reflective of the tissue bulk under investigation. Depending on how parameters are discretised, MRF dictionaries can become very large and the computational burden of matching the signal to the dictionary expands. Irrespective of the approach taken to estimate model parameters, the parameters are depicted as spatially resolved maps (i.e., model is applied at each pixel/voxel location), which are quantitative images.
MRF signals have mostly been generated using the classical Bloch equations and other approaches not involving fractional calculus [18,19]. Recently, the use of the time-fractional Bloch equations was investigated [20]. The MRF approach provides a platform for exchanging models with an associated adjustment of the parameter space. Building on previous work on anomalous relaxation in MRI [11], our aim here was to investigate further the link between the time-fractional order parameters and expected variation in tissue microstructure. The overall motivation being that accurate delineation of microstructurally distinct regions of the human cerebral cortex in individuals, i.e., known as cortical parcellation, is fundamental for understanding structure-function relationships in the brain [21]. Additionally, specific variations in model parameters may shed light on aging and how different diseases affect cortical regions [22,23], and improve delineation of abnormal tissue in the surgical setting [24,25]. Current methods in the cortical parcellation area mostly involve the assessment of changes in MRI relaxation times (such as 1 and 2 ) [26][27][28][29], and it is becoming increasingly evident that these model parameters are unlikely capable of distinguishing interareal structural variations throughout the whole human cerebral cortex [30][31][32][33]. In what follows we outline our approach of using the time-fractional Bloch equation in the MRF context for parcellating the cerebral cortex in individuals.

Bloch equations
In this subsection we provide the formulation used in the study based on a previously published model [11]. We first describe the general, time-fractional order Bloch model, which can be simplified to the classical case. The equations describing the evolution of magnetisation within a spin system and then related to an MRI signal.

Time-fractional order model
We consider the case of spin-spin relaxation in an inhomogeneous magnetic field due to microscale variations in tissue microstructure. In the rotating frame a residual frequency offset presents [11], Δ = 2 Δ where Δ is in units of Hz, and 2 is replaced by 2 * and strictly, relaxation times satisfy 1 ≥ 2 ≥ 2 * and measured in units of s. Whilst different approaches may be considered to arrive at a time-fractional order Bloch equation [15], here the approach of fractionalising the time derivative was taken: where constants 1 and 2 are incorporated to preserve units, 0 in units of A/m is the net magnetisation produced by the MRI scanner. We should note that the z-component magnetisation, where Γ denotes the Gamma function. The solution to (1), as described in [11], takes the following form: is the two-parameter Mittag-Leffler function, and by definition ,1 ( ) ≡ ( ) and 1 ( ) = . Noticeably, ( ) and ( ) involve a different time-fractional order than ( ). This is because from a physics perspective, different processes drive magnetisation change for these components (i.e., 1 versus 2 * effects). Since 1 ≫ 2 * , the time scale over which ( ) recovers to 0 is also quite different to the time scale over which ( ) and ( ) tend to zero. This becomes important later in the context of MRI data collection.
The condition 1 ≥ 2 ≥ 2 * holds in the case of the integer order model. In the time fractional case, we can show that It is possible to additionally set Δ = 0 in (3), which would then be the case of spin-spin relaxation in the absence of field inhomogeneities. Note, according to the Larmor equation (i.e., Δ = 42.578 × 10 6 Δ for hydrogen spins where Δ is the induced change in magnetic flux density), an induced change in field results in an induced shift in frequency. During MRI data acquisition it is possible to apply a so-called radio frequency refocusing pulse which ensures signals are not de-phased at the time of signal acquisition. With the use of such an approach 2 * in (4) would have to be exchanged by 2 and Δ = 0 can be assumed. The reduced form time-fractional Bloch equations has been used previously as well [35,36].

From magnetisation to an MRI signal
The spin system considered in MRI results in a net magnetisation vector which can be considered to precesses in each MRI voxel, known as a volume pixel. Since the reference magnetisation, 0 , is defined in the z-coordinate direction, precession of spins occurs in the plane perpendicular to the zdirection. That is, in the reference frame precessional freqeuncy exists for ( ) and ( ), and not for ( ) (see in (3) that a frequency shift is only present on ( ) and ( ), and for additional informatio refer to [11]). This means the signal detection system, which is sensitive to a specific frequency band, produces a signal due to ( ) and ( ) alone. Traditionally, the induced signal in the induction coils used to collect the MRI signal are converted to a magnitude signal, such that: which after the subsitution of (3) conveniently becomes: where 0 = √ 2 (0) + 2 (0) is the signal when = 0. After careful inspection, we find that (6) is insensitive to ( ), suggesting the z-component of magentisation behaves independently of ( ) and ( ). However, this is not stricly true. In an MRI system, magnetisation at = 0 is preserved, meaning 0 = √ 2 (0) + 2 (0) + 2 (0). Rearranging of this relationship leads to two definitions for 0 , namely as in (6) and additionally 0 = √ 0 2 − 2 (0). In the MRF context, the 0 term involving 0 and (0) ensures estimation of all models parameters can be approximated by matching (6) with the experimentally acquired MRI signal. We should point out that in the 2 * regime the temporal magnetisation vector is complex valued, implying that a phase is also generated during MRI data acquisition, which takes the value = × based on (4). This phase can be interpreted as a shift in frequency arising due to a microscale variations in tissue magnetic properties inducing small magnetic field variations. Note, according to the Larmor equation, = 2 × 42.578 × 10 6 × , phase increases with and . We should additionally point out that ( ) is real valued, and it is not influenced by . Hence, one may expect static field inhomogeneities introduced by the sample within MRI voxels to only influence ( ) and ( ), and not ( ).

Parameter estimation using MRF
Suppose the MRI instrument collects signals, ( ) . Then one approach of generating model parameters would be to fit the signal model, described by (3) and (5), to the acquired MRI temporal signal, ( ). This approach requires a non-linear fitting algorithm with good initial parameter guesses to be able to achieve reasonable parameter estimates, especially when MRI data contains noise, as is the case normally [11].
In MRF a dictionary matching approach to parameter estimation is considered. Suppose the parameter space is known, and it can be represented using discrete values. The MRF dictionary is often represented as a matrix, where each column defines a specific parameter, and each row is a unique set of parameters. The key to each dictionary entry is the unique MRF signal, ( ), which can then be simulated using (3) and (5) for each parameter combination [17]. The best matched ( ) to ( ) is considered to have the parameters which best parameterise the voxel. After performing the matching for each voxel, where commonly the metric to match signals is the dot product between the time series signals, i.e.
( ) ⋅ ( ) over parameters P, the resultant parameters are depicted as spatially resolved maps, sometimes referred to as quantitative images. Based on (3), each voxel can be parameterised in terms of 1 , 2 * and . Time fractional exponents and are assumed to lead to additional parameteric maps containing information relating to tissue microstructure and constituents. For convenience 1 and 2 are assumed to equal 1.
The first step in generating the MRF dictionary is to decide on practical bounds for each of the parameters, often constrained by existing knowledge. For example, it is known that gray and white matter 1 values in the brain range between 500ms to 3000ms, and similarly bounds can be placed on 2 * and . The fractional exponents by definition have to satisfy , ≤ 1, and typically they are close to 1. Here, we opted to create the MRF dictionary using the following parameter ranges: {500 One may consider the MRF signal matching approach as a constrained parameter estimation problem. The MRF dictionary is essentially a discrete representation of the parameter space defined through certain combinations of individual parameters. The number of parameter sets in the dictionary can become very large rapidly, as governed by the increments between parameters. Also, the creation of the discrete parameter space from which signals are simualted is a combinatorial problem, wherein physical limits can be adopted to reduce MRF dictionary size.

Relating MRI data to the Bloch model
The flip angle ( for simplicity) used in the MRI data acquisition equation 'flips' the (0, 0, 0 ) magnetisation vector by an angle with respect to the z-axis. This, when applied at = 0, results in a new z-component magnetisation, such that (0) = 0 . Substitution of this value into 0 = √ 0 2 − 2 (0) yields 0 = 0 , which is the term governing the MRI signal amplitude in (6). After the application of an initial flip of the magnetisation vector, the signal is collected at = . However, the next application of the flip angle occurs at = , at which point, according to (4), left over zcomponent magnetisation is ( ). This value becomes the (0) for the next TR cycle, or MRF repetition. In this way information on 1 is encoded into 0 . Note, changes in flip angle lead to a nonlinear modulation of the MRI signal, essentially bringing more complexity into the MRF signal repetitions. The matching of the simulated and acquired signals (as a function of MRF repetitions over TRs) leads to a signal which best matches the measured MRI signal. The parameters used to generate that simulated signal are then used to create parametric maps in a pixel-by-pixel manner.

MRI data collection
Six mixed gender volunteers between 27 and 35 years of age participated in the study. The MRI sequency was implemented by SM and used to acquire data using the 7T whole-body MRI research scanner (Siemens Healthcare, Erlangen, Germany) located at the Centre for Advanced Imaging, University of Queensland, Brisbane, Australia. Data acquisition involved a 3D echo planar imaging (EPI) MRF sequence [37] with 1000 frames made up of 3D MRF images with the following parameters: repetition time (TR) = 41-99 ms, flip angle (FA) = 10-77°, echo time (TE) = 12-48 ms, partial Fourier phase = 6/8, voxel size = 1.4 × 1.4 × 1.4 mm, and matrix size = 142 × 142 × 88. We used GRAPPA parallel imaging [38] in both phase encoding (with acceleration factor = 3 and reference lines = 36) and slice encoding directions (with acceleration factor = 2 and reference lines = 12). Chemical-shift-selective (CHESS) fat saturation technique [39] was used to reduce common artefacts observed in EPI sequences at high field scanners [40]. The sinusoidal FA pattern used for the MRF acquisitions, shown in Figure 1, was assumed from [41]. Abrupt FA changes at the final MRF frames were added for increased sensitivity of the MRF signals to transmit field variations, as suggested previously [42].   We additionally used pseudo-randomised patterns of TE variation suggested by Rieger et al. [43], see Figure 2, to improve signal-to-noise ratio as a result of the large number of small TE values. The use of small TE values also led to a reduction in total acquisition time, as TR (the sum of which scales with the total data acquisition time) can be made proportionally small with resepct to TE. The alternating TE pattern has also been shown to increase sensitivity of MRF signals to 1 and 2 * variations [43], the positive outcome of which is improved estimation of these values.

Model selection and estimation error
Traditionally the integer order Bloch equations (i.e., (4)) have been used to generate simulated MRF signals based on predefined sets of discrete parameters. By replacing the classical Bloch equations with the time-fractional counterpart (i.e., (3)), the number of model parameters increases by three, namely , and , assuming 1 = 2 = 1. According to the Akaike information criteria [44], for example, the residual of the fit would have to be penalised due to an increase in the number of model parameters. As such, the residual error of fit using the time-fractional model is expected to be smaller than using the classical approach, and a penalty associated with increased numbe of model parameters has to be considered. Nonetheless, we are particularly interested not only if the time-fractional Bloch model is able to better fit the data, but also if the additional model parameters contain new information. The ability to provide information on cyto-architecturally different brain regions based the timefractional Bloch equations could suggest that time fractional order plays a key role in MRI signal formation.

Results
We provide MRF simulation results based on the time-fractional Bloch equations, followed by results on the separability of different cortical regions in the human brain based on time fraction Bloch model parameters.

Expected changes in the MRF signal
The MRF signal follows a pseudo-random pattern as defined by the acquisition protocol flip angle (i.e., ), TE and TR. The plots provided in Figure 2 show how the MRF signal based on the acquisition protocol parameters provided in Figure 1 evolve over MRF acquisition repetitions. The signal is acquired at = TE, and ( ) is allowed to evolve until t = TR, and then the sequence is repeated. The different plots consider changes in 1 , 2 * , , and , and in each case a freqeuncy shift of = 10Hz. Distinct changes in 1 (Figure 2(e)) and 2 * (Figure 2(f)) lead to specific changes in the MRF signal over repetitions. We may also notice that a change in (Figure 2(b)) leads to a change in ( = ), essentially modulating the amount of magnetisation available for each repetition. The parameter acts on the transverse magnetisation, and thus on ( ), the transverse magnetisation derived MRF signal (Figure 2(c)). Changes in both 2 * and appear to amplitude modulate ( = ).

Time-fractional Bloch model parameter sensitivity
Time fractional Bloch equations simulations were performed based on the parameter ranges provided in Section 2.2. For each simultion, i.e. a data point, a random value for 1 , 2 * , , and were chosen from their respetive intervals. Additionally, a random change was applied to each parameter and assessed in a one-by-one manner, and for the sake of this analysis, repeated 500 times (i.e., 500 data points were generated in each case). The error applied to each parameter was limited to a maximum of 25% with respect to the actual parameter value. This process led to the generation of parameter sensitivity maps in view of the angle produced by the dot product matching, as shown in Figure 3. It appears from the results provided that the time-fractional Bloch model, which includes 1 , 2 * , , and , overfits the MRF signal, see wide ranging values for 2 * , and . Values for 1 and appear not to be overfitted, and good sensitivity to these parameters can be achieved. See slope on the fit and also the clustering of data points about the linear regression line.
Noting that 2 * is a physical parameter, and findings from Figure 3 suggest model degeneracy due to an interplay between and , as indicative of results shown in Figure 2 wherein both of these parameters amplitude modulate ( = ). To overcome the overfitting problem, we may consider two options; set = 1 or = 0. The former choice causes a fundamental problem, since subsitution of = 1 into (6) leads to a signal insensitive to changes in . As such, we opted to investigate the option of settting = 0 and leaving as a free model parameter. In Figure 4 parameter sensitivity results when = 0 are provided. The results suggest that parameter insensitivities concluded from Figure 3 can be overcome by setting = 0. Judging from the regression line slope, the best sensitivity is achieved for , followed by and 1 , and lastly 2 * . These results are in view of the paradigm described in Figure 1, and may not be generalisable across different MRF acquisition protocols wherein differnet TEs and TRs may be set.

Parameter selectivity to different cortical regions in the human brain
For the purpose of being able to appreciate how data were acquired, Figure 5 illustrates an MRF image slice in the human brain and MRI signal evolution over MRF repetitions. At each location in the brain in the image, the MRF matching based on the dot product metric was applied and parameters associated with the best matched signal were deemed to reflect location specific tissue parameters (i.e., 1 , 2 * , and ). Figures 6 and 7 illustrate the spatially resolved maps for each of the time-fractional Bloch equations paramters for an example slice at different locations for two of the participants. Figure 8 provides the boxplot results for matching of the MRF simulated signal with the acquired MRI signal based on the interger order Bloch equations and time fractional order counterpart. It has previously been demonstrated that 1 and 2 * are not selective to cortical regions in a way full cortical parcellation of the human cerebral cortex can be achieved in individuals. We also found that 1 and 2 * mapped using the MRF protocol is unable to differentiate the cortical regions identified in Section 2.6. However, we did find both mapping of and provides specificity to different cortical regions. Figure  9 are violin plots for and over the various regions categorised into four groupings: somatosensory cortical region (BA1, BA2, BA3a and BA3b), the primary motor (BA4a and BA4p) and pre-motor (BA6) cortical region, visual cortex region (BA17 and BA18) and two adjacent Broca cortical regions (BA44 and BA45). Based on the notch plot, was found to be significantly different in the somatosensory cortical region and the primary and pre-motor cortical region. Interestingly, was significantly different in the sub-regions of the four cortical areas studied. This finding is very interesting, as differentiation of the Broca regions are usually among the most challenging parcellation tasks in neuroscience.
In Figure 10 the relationships between time fractional Bloch equations parameters are investigated in the entire human brain. We perfomed tests when was a fitted parameter. In this case, see Figure  10(a), and values take on wide ranging values, and appears to be linked with , Figure 10(b). By setting = 0, a vastly different relationship between and presents (compare Figures 10(a) and 10(c)). A systematic trend between and when = 0 cannot be discerned.

Discussion
The aim of the research was to establish the utility of the time-fractional Bloch equations in magnetic resonance fingerprinting for brain imaging studies. The classical, integer order in time Bloch equations, are known to produce a substantial residual in the MRF fit, and residuals contain information not explained by classical model parameters [52]. The incorporation of time-fractional exponents led to a reduction in the fitting error (see Figure 8). However, a dependency between fractional exponent and Δ was found, resulting in overfitting (see Figures 2 and 3). When Δ = 0 in the time fractional Bloch equations, overfitting does not occur and exquisite spatially resolved maps of and can be obtained (see Figures 6 and 7). This approach allowed us to separate out different cortical regions in the brain, with providing better results than (see Figure 9).

MRF parameter discretisation and matching
Another study considered time-fractional order Bloch equations in an MRF implementation [20]. The method was applied in a 1 / 2 phantom and using the time-fractional model better estimates of 1 and 2 were generated in comparison with the integer order Bloch equations implementation. Notably, the authors considered the case of 2 relaxation, i.e., no signal dephasing produced based on the spinecho acquisition protocol. The parameters were discretised using the following choices: 1 values in the range [100, 4500] in steps of 100ms; 2 values in the range [10,1000] in steps of 10ms; both and were in the range [0.96, 1.1] in steps of 0.01. The number of MRF repetitions was 600. It was found that 1 and 2 estimates improved when and were free parameters, and it was additionally reported that > .
(a) (b) Figure 8. A comparison between the matching using the integer order and the time fractional order Bloch models. Shown are the (a) root-mean-squared error (RMSE) and (b) the angle (i.e., in degrees) produced by the dot product for matching the simulated MRF signal with the one acquired using the MRI instrument. The box plots conclude that the fractional order Bloch model was able to produce significantly smaller RMSE and (p < 0.05).

SOMETOSENSORY PRIMARY & PRE-MOTOR VISUAL BROCA
(a) (b) Figure 9. Violin plots for and including notches at the (p < 0.05) significance level across the 11 cortical regions considered. The headings identify different cortical areas consisting of different regions, for example adjacent cortical areas BA17 and BA18 in the visual cortex have been evaluated. The time fractional Bloch model parameter (a) was not able to discern regions (p < 0.05) in the visual and Broca cortical regions, whereas (b) had significant differences between sub-regions in the four cortical regions evaluated.
The MRF dictionary contains a set of discrete parameters, and the finer the increments the larger the time required to obtain a match. As such, there is a trade-off between dictionary size and computational resource requirements to make a signal match. We chose to have a coarse-grained approach: 1 and 2 * increments of around 1.5-4%, whereas , and Δ as much as 10% increments. Results presented in Figure 4 suggest that for a 10% change in 1 one can expect around 1.6% change in matched angle (assuming 90 is the maximum). For 2 * at the same level of change we expect about half as much change in the dot product produced angle between MRF simulation and MRI data acquisition over repetitions. Interestingly, has the largest influence on the angle, which is about double that of 1 , and is similar to 1 . Based on these findings, in the future consideration should be made on how the different model parameters are discretised. Ideally, the steps chosen for each parameter should approximately produce the same change in the angle between the MRF repetitionbased acquired and simulated signals.

Role of and
In the classical case = = 1, which reduce the time-fractional Bloch equations to their integer order form. We opted to perform the study using distinct values for and , as it has been shown that these values are unlikely to be the same [20,35,36]. These previous studies have suggested > , and additionally it was proposed that = 1, resulting in mono-exponential recovery of magnetisation to 0 [35,36]. Interestingly, our brain imaging findings suggest that > does not hold in the brain, refer to Figures 6 and 7, and also should not be set to 1. Out of interest, we provided sensitivity results for when = 1, see Figure 3(f), and found this choice not to provide a benefit for model fitting.
It is well established that both 1 and 2 have a level of frequency reliance, which depends on the underlying composition of the material imaged, and likely scales as a power law [53]. This effect is particularly pronounced in the presence of lipid/protein structures, where dipole-dipole coupling, for example, frequents [54]. The brain consists of large proportions of proteins and lipids, myelin being an example of highly concentrated and organised lipid structures in the brain. It is also known that lipids generally have magnetic properties which influence the MRI magnetic field responsible for the net magnetisation, 0 . Our human brain results in Figure 10 showed a potential relationship between and Δ , the latter being the field change induced frequency shift in the MRI signal. We found that the use of both parameters led to overfitting (Figure 3 versus Figure 4), suggesting one of these two parameters is sufficient in explaining the trend in the MRI signal. Simulation findings presented in Figures 3 and 4, and subsequently in the human brain, see Figure 10, imply that and Δ vary together. In view of previous findings on 1 and 2 frequency dependence and noting that the difference between 2 and 2 * is attributed to tissue induced magnetic field inhomogeneities, our results, interestingly, suggest that captures information on induced field change. That is, an increase in Δ corresponds with a decrease in , particularly above Δ = 20 , and remains to be further evaluated for small Δ values. Ideally an analytic expression linking the two would provide the best insight into the relationship between these two parameters. Figure 10 results additionally suggest that = is unlikely, and setting = 1 does not seem appropriate for brain studies. Whilst in this work we have not attempted to demonstrate a clear relationship between 2 and 2 * , it would be interesting to find 2 * = 2 . This is purely an observation based on 2 and 2 * values reported for the brain, and future work would need to investigate carefully whether such a relationship is mathematically plausible.

Cortical parcellation
Parcellation of the human cerebral cortex in individuals is challenging and a robust non-invasive imaging method of achieving this goal has not been proposed to date. The time-fractional Bloch equations appears to provide good sensitivity to different cortical regions through model parameters, as illustrated in Figure 8. In our study the time-fractional exponent of the transverse components of the magnetic field, , led to better cortical area differentiation than , the time-fractional exponent influencing longitudinal magnetisation evolution. Previously it was shown that the integer order Bloch equations within an MRF framework result in signal residuals depicting trends useful in delineating cortical regions of distinct cyto-and myelo-architectures [52]. The hypothesis that the cyto-architectures of cortical regions lead to induced magnetic field changes distinct to regions appears valid, based on the specificity of , a plausible proxy for cortical tissue induced field changes via Δ .

Conclusions
Mathematical model-based approaches of extracting information from MRI data provide an important avenue for producing quantitative, parametric maps, reflecting tissue properties. Integer order models have found excellent utility in MRI-based parametric mapping, and fractional calculus, whilst explored, has had limited impact to date. Our time-fractional order Bloch equations magnetic resonance fingerprinting results suggest that non-integer order approaches can lead to a better explanation of the trends in the magnetic resonance fingerprinting signal. Moreover, the time-fractional exponents, one associated with magnetisation recovery and the other with transverse magnetisation loss, can provide new insights in human brain studies involving tissue specific parameter estimations. Our example application involved the parcellation of the human brain using time-fractional exponents of the Bloch equations.