# A Statistical Robust Approach to Design Parallel Transmit Radiofrequency Excitations in MRI

Authors: Vincent Gras; Franck Mauconduit; Nicolas Boulant
Journal: Concepts in Magnetic Resonance Part A
 (2020)

---

## Abstract

In MRI, at ultrahigh field, the use of parallel transmit radiofrequency (RF) arrays is very beneficial to better control spin excitation spatially. In that framework, the so-called “universal pulse” technique, proposed recently for head imaging at 7 tesla, gives access to “plug-and-play” nonadiabatic solutions exhibiting good robustness against intersubject variations in the resonant transmit fields. This new type of solution has been defined so far as the result of numerical pulse optimizations performed across a collection of RF field maps acquired on a small population sample (pulse design database). In this work, we investigate an alternative universal pulse design approach in the linear small tip angle regime whereby the database of RF field maps is first transformed into a second-order statistical model and which then exploits a statistical robust design formalism for the optimization of the RF and magnetic field gradient waveforms. Experimental validation with an eightfold transmit RF coil for 7 tesla brain imaging shows that this new approach brings some benefit in terms of computational efficiency. Hence, for a design database composed of 35 maps, the computation time initially of 50 min could be reduced down to 3 min. The proposed statistical approach thus enables integration of large databases, presumably necessary to ensure robust solutions. Finally, it provides means to compute flip angle statistics and, along with it, simple performance metrics for quality assurance (RF pulse performance) or guidance in the optimization of TX array architectures.


---

## Body

## 1. Introduction

In magnetic resonance imaging at high magnetic field, the apparition of standing wave effects [1–3] prevents uniform excitation of the spins across extended anatomical regions. For brain imaging, these effects are clearly visible at 3 tesla (T) and are often deemed unacceptable at field strength >4 T. Parallel transmission (pTX) is a promising technology to mitigate this problem as it provides degrees of freedom to “shim” the RF field distribution [4] (RF shimming). The transmit SENSE technique [5–7], enabled by the pTX technology, furthermore allows reducing the duration of multidimensional pulses [8] while preserving spatial definition. It thus offers one more solution to tackle the transmit field heterogeneity problem. The rationale here is to regard the problem of achieving a uniform magnetization profile as one particular problem that multidimensional RF pulses can solve, whereby the target magnetization profile is spatially uniform. The design of such pulses, here referred to as flip angle (FA) shimming, plays now an important role at UHF magnetic resonance to enable uniform spin excitations. Thus, the so-called fast-kz or spokes pulses [9, 10] and the kT-point pulses [11] are well-known slice-selective and nonselective FA shimming strategies taking advantage of the parallel transmission technology to produce highly uniform excitation profiles at a reasonable cost in terms of pulse duration as compared to RF shimming.However, in addition to being nonuniform in space, the transmit RF field (B1+) also depends on the geometry and electric properties (dielectric permittivity and conductivity) of the sample and thus can vary significantly across individuals. Thus far, the usual technique to cope with this is to characterize the actual B1+ profiles using MRI-based measurements and to tailor RF shimming to this particular TX field distribution. This obviously can be time-consuming, in particular when the number of TX channels is high [12]. One alternative approach to handle the variability of the TX field distribution, called universal pulse (UP) [13], is to replace the subject-specific TX field distribution with a collection of such TX field distributions acquired on a database of subjects. The rationale here is to achieve, simultaneously across all subjects of the database, good excitation performances, and assume that this property generalizes to the population that the database is supposed to cover (e.g., the adult population). The validity of this approach has been demonstrated experimentally for various applications and hardware setups [14–16]. In those studies, a database of 10 to 20 individuals was deemed sufficient to yield good performances. However, the amount of data remained too limited to analyze quantitatively how pulse performance changes with the design database, and in particular with its size. Nevertheless, it is reasonable to assume that UP performance increases with the size of the design database; conversely, if a too small database is used to compute UPs, there is a risk of “data overfitting,” i.e., that the solution displays a higher FA shimming performance on the design database than its actual performance on the entire population.To better control the actual performance of UP on the entire population, we propose in this work to model the transmit RF field as a stochastic process, whereby every new subject (giving rise to a newB1+ map) is a realization of this process. This representation appears very useful as it allows formulating the action of a FA shimming pulse in statistical terms. From the RF pulse design perspective, it also provides an adequate framework to apply the so-called statistical robust design theory, proposed by Alotto et al. [17], applied to the design of electromagnetic devices. For the latter application indeed, manufacturing errors of the components or uncontrollable environmental conditions such as temperature or humidity can lead to uncertainties in the response of the system. Applying classical optimization to such a problem can severely degrade the actual performance of the solution. A robust design approach mitigates this problem by accounting for the randomness of the response of the system in the optimization process using an appropriate statistical model. This idea can naturally be extended to the RF pulse engineering problem.We thus present in this work a statistical representation of the TX field distributionB1+ in which, for every position r, B1+r is assimilated to a Gaussian random process. We then derive a method to apply the aforementioned statistical robust design method to the design of FA shimming UPs, which uses the second-order statistics of B1+r instead of the collection of B1+ maps as input parameter. An experimental validation of the method is provided with data acquired with an eightfold transmit RF coil for 7 tesla head MRI. We demonstrate in particular (1) the utility of statistical robust design approach to boost computational efficiency, (2) the capacity of the method to enable the integration of a large database TX field maps, and (3) the utility of a statistical framework to compute flip angle statistics and to derive simple performance metrics for quality assurance.


## 2. Theory

### 2.1. Statistical Modelling of the Resonance TX Field

LetB1,c+r (tesla/volt) denote the spatially dependent (r referring to the position) resonant transmit RF magnetic field of the cth resonator of the Nc-fold transmit array and γ=42.57×106 Hz/T the gyromagnetic ratio of the proton. The so-called control field components ω1,j=2πγB1,j+ (rad/s), 1≤j≤Nc, may then be arranged in a vector ω1r=ω1,1,…,ω1,Nct∈ℂNc×1 called the control field vector.The control field vector being subject dependent is viewed as a stochastic variable where each new subject gives rise to a realization ofω1. For every location, we assume that the mean value ω1r and the covariance matrix (∗ for the Hermitian conjugate) ω1r−ω1rω1r−ω1r∗ of the control field vector is well defined and we denote them, respectively, by μr∈ℂNc×1 and Cr∈ℂNc×Nc. In the following, the control field ω1r is assumed to have Gaussian statistics, i.e., ω1r∼Nμr,Cr.Letω1N,…,ω1N be N independent control field distributions (in practice TX field maps measured in N different subjects). If N is sufficiently large, this dataset can be exploited to estimate the mean μr and covariance matrix Cr of the underlying random process by computing the sample mean and covariance matrix for this ensemble:(1)μr≃1N∑n=1Nω1nr=μ^v,Cr≃1N−1∑n=1Nω1nr−μ^rω1nr−μ^r∗.For the covariance matrix, we shall use the Ledoit–Wolf [18] correction of the sample covariance estimator when the condition N≫Nc2 (Nc2=64 in this study) is not satisfied.Given a statistical model forB1+, we now show that this information is sufficient to compute universal FA shimming pulses and to study the robustness of such pulses with respect to the intersubject variability of B1+. We first propose a formalism to describe the action of FA shimming pulses. Although a generalization would be possible for large flip angles, this theory for now is limited to the small tip angle regime where there exists a linear relationship between the RF excitation and the spin’s tip angle [19].


#### 2.1.1. Matrix Representation of the Spatial Distribution the Transmit RF Field

As the control field distribution is typically obtained by MR acquisition, it is common in practice thatω1 is not defined at each voxel of the underlying MR acquisition (existence of regions where the reconstruction of B1+ is not possible). We thus define by ℛ the set of voxel positions r where a given realization of ω1r is known (ℛ is subject and session dependent). In the following, we shall use the following matrix representation of the spatial distribution of the control field:(2)Ω1ℛ=ω1rr∈ℛ∈ℂNc×M,where M=ℛ is the number of positions contained in ℛ. When there is no ambiguity about the choice of the region of interest, we shall simply use the notation Ω1.


### 2.2. Flip Angle Shimming Pulse

Letut∈ℂ1×Nc and kt∈ℝ1×3 be the RF and k-space trajectory waveforms of a given FA shimming pulse. Let T be the duration of the RF and magnetic field gradient (MFG) waveforms and ω0 (rad/s) be the resonance offset at location r due to the static field nonuniformity. Let us use below the notation Pr,t=uteikt⋅r. In the small tip angle approximation, the effective (time-averaged) control field created by this RF and magnetic field gradient (MFG) pulse pair is given by [8](3)ω¯1P,r=P¯rω1r,where(4)P¯r=1T∫t=0Tuteikt⋅reiω0rT−tdt.Note that the flip angle created by the pulseP is simply ϕP,r=Tω¯1P,r. As a result, the mean, μω¯1, and variance, σω¯12, of the effective control field ω¯1 are given by(5)μω¯1r=P¯rμr,σω¯1r2=P¯rCrP¯r∗.In this model, the effect of off-resonances is taken into account but is not modelled as a random processω0r. Equation (4) hence can be implemented as such if the static field offset is mapped in a subject-specific manner. Alternatively, one can assume ω0=0, i.e., that the static field offset has a negligible influence on pulse performance, which is acceptable for short (Tω0≪1) pulses.In the following, we present a method to compute universal FA shimming pulses which is entirely based on a second-order statistical model for the control field. Along with this statistically driven framework, we provide a method to evaluate the performance of a universal FA shimming pulse.


### 2.3. Statistical Robust Design of Universal Pulses

LetP be a candidate FA shimming pulse to create a unitary effective control field ω¯1. Let us denote by(6)δP,r=ω¯1P,r−1,the discrepancy between the unit target and ω¯1P,r, and by ΠPΔ,r the probability that δP,r exceeds a certain error level Δ>0. To optimize P, we can exploit here a statistical robust design (SRD) approach [20]. Let us introduce a new design parameter 0≤θ≤1, called SRD tolerance, and define the error level ΔP,θ as(7)ΔP,θ=minΔ≥0;ΠPΔ,r≤θ∀r∈ℛdesign,where ℛdesign denotes the region of interest for the design of the pulse (any position outside r being therefore ignored in this metric). The SRD problem may then be posed as the minimization problem:(8)minPΔP,θ.Denoting now byFδP,r the cumulative distribution function of δP,r, the condition ΠPΔ,r≤θ gives(9)FδP,rΔ≥θ.Hence, the problem described by equation (8) takes the equivalent form:(10)minPmaxr∈ℛdesignFδP,r−1θ,where FδP,r−1 denotes the inverse function of FδP,r. A useful simplification of this optimization problem is then obtained by replacing the maximum (ℒ∞ norm) by the ℒq norm for a certain choice of q≥1. This leads to the following new pulse optimization problem:(11)minP∑r∈ℛdesignFδP,r−1θq.If, now, for a given pulse candidateP, the following two conditions are satisfied for all positions r in ℛdesign:(12)varω¯1P,r≪ω¯1P,r2,ω¯1P,r−12≪varω¯1P,r,then we can show (see Appendix A) that the pulse optimization problem (11) can be further simplified to take the form:(13)minP∑r∈ℛdesignω¯1P,r−12q/2,where the dependence on the second-order B1+ statistics can be made explicit:(14)ω¯1P,r−12=12P¯rCP¯r∗+P¯rμr2−12+12P¯rCP¯r∗.Retrospectively, the link with the database-driven objective function (simultaneous design across all subjects of the database ofB1+ maps), denoted here by fUP, can be explained as follows. The latter can be defined as the mean across the design database of the mean squares error of the effective control field:(15)fUPP=∑n=1Ns∑r∈ℛdesignP¯rω1nr−12.Now, if we assimilate sample mean and statistical mean, we obtain(16)fUPP=∑r∈ℛdesignP¯rω1r−12=∑r∈ℛdesignω¯1P,r−12,which leads us to the same expression as the one obtained in equation (13) for q=2. This shows the near equivalence (note however that the objective function defined in reference [15] is proportional to the average (ℒ1 norm) of the FA-NRMSE rather than the ℒ2 norm and that the brain region for the evaluation of the FA-NRMSE was subject-specific rather than fixed as in equation (15)) of the database-driven and the statistical approaches when sample mean and statistical mean can be assimilated. This is acceptable in particular when the design database is large enough to be representative of the whole targeted population.


### 2.4. Derivation of a Lower Bound for the Coefficient of Variation of the Flip Angle

Letαr be a solution to the problem of minimizing the variance of the control field at position r among the RF shim solutions that return on average, up to a phase term, a unitary control field:(17)αr=argmina∈ℂ1×Ncaω1r−12subject toaω1r=1.In AppendixB, an explicit description of αr is obtained under the condition that(18)βr≜μrC−1rμr∗≫1.The expression forαr then reads(19)α=C−1μ∗β,and the variance of the value of the objective at the optimum is given by(20)αrCrαr∗=1β.Furthermore, it is shown in AppendixC that the locally optimal solution αr sets a performance limit for the coefficient of variation (CV) of the effective control field:(21)gr=σω¯1rμω¯1r,which reads(22)gr2≥αrCrαr∗.Hence, under the conditionβ≫1, the lower bound for the CV of the effective control field, which we may also call ultimate FA-CV in the linear STA regime, and below denoted by g∞, reads(23)g≥1β≜g∞.The conditionβ≫1, on the other hand, has a clear physical interpretation: the mean control field magnitude μ dominates its standard deviation. The validity of this condition in practice will be tested subsequently using experimental data.It is also of interest to define the ultimate FA-CV—in the sense of equation (23)—for the circularly polarized (CP) mode of the TX array. In this mode, a fixed relationship between the RF signals applied to the TX array ports is imposed in order to reproduce the behavior of a single port volume coil. The CP mode can be defined with the vector a⌣∈ℂ1×Nc which describes the way the channels of the TX array interfere. The mean and variance of the CP mode control field, μ⌣,C⌣, is then given by μ⌣=a⌣μ and C⌣=a⌣Ca⌣∗. Application of equation (23) on the CP mode control field statistics gives(24)g⌣∞=1β⌣,where(25)β⌣=a⌣μ2a⌣Ca⌣∗.In what follows, we present an application of the proposed statistical model of the intersubject variability of the control field and their implication in terms of pulse design and flip angle shimming performance assessment.


## 2.1. Statistical Modelling of the Resonance TX Field

LetB1,c+r (tesla/volt) denote the spatially dependent (r referring to the position) resonant transmit RF magnetic field of the cth resonator of the Nc-fold transmit array and γ=42.57×106 Hz/T the gyromagnetic ratio of the proton. The so-called control field components ω1,j=2πγB1,j+ (rad/s), 1≤j≤Nc, may then be arranged in a vector ω1r=ω1,1,…,ω1,Nct∈ℂNc×1 called the control field vector.The control field vector being subject dependent is viewed as a stochastic variable where each new subject gives rise to a realization ofω1. For every location, we assume that the mean value ω1r and the covariance matrix (∗ for the Hermitian conjugate) ω1r−ω1rω1r−ω1r∗ of the control field vector is well defined and we denote them, respectively, by μr∈ℂNc×1 and Cr∈ℂNc×Nc. In the following, the control field ω1r is assumed to have Gaussian statistics, i.e., ω1r∼Nμr,Cr.Letω1N,…,ω1N be N independent control field distributions (in practice TX field maps measured in N different subjects). If N is sufficiently large, this dataset can be exploited to estimate the mean μr and covariance matrix Cr of the underlying random process by computing the sample mean and covariance matrix for this ensemble:(1)μr≃1N∑n=1Nω1nr=μ^v,Cr≃1N−1∑n=1Nω1nr−μ^rω1nr−μ^r∗.For the covariance matrix, we shall use the Ledoit–Wolf [18] correction of the sample covariance estimator when the condition N≫Nc2 (Nc2=64 in this study) is not satisfied.Given a statistical model forB1+, we now show that this information is sufficient to compute universal FA shimming pulses and to study the robustness of such pulses with respect to the intersubject variability of B1+. We first propose a formalism to describe the action of FA shimming pulses. Although a generalization would be possible for large flip angles, this theory for now is limited to the small tip angle regime where there exists a linear relationship between the RF excitation and the spin’s tip angle [19].


### 2.1.1. Matrix Representation of the Spatial Distribution the Transmit RF Field

As the control field distribution is typically obtained by MR acquisition, it is common in practice thatω1 is not defined at each voxel of the underlying MR acquisition (existence of regions where the reconstruction of B1+ is not possible). We thus define by ℛ the set of voxel positions r where a given realization of ω1r is known (ℛ is subject and session dependent). In the following, we shall use the following matrix representation of the spatial distribution of the control field:(2)Ω1ℛ=ω1rr∈ℛ∈ℂNc×M,where M=ℛ is the number of positions contained in ℛ. When there is no ambiguity about the choice of the region of interest, we shall simply use the notation Ω1.


## 2.1.1. Matrix Representation of the Spatial Distribution the Transmit RF Field

As the control field distribution is typically obtained by MR acquisition, it is common in practice thatω1 is not defined at each voxel of the underlying MR acquisition (existence of regions where the reconstruction of B1+ is not possible). We thus define by ℛ the set of voxel positions r where a given realization of ω1r is known (ℛ is subject and session dependent). In the following, we shall use the following matrix representation of the spatial distribution of the control field:(2)Ω1ℛ=ω1rr∈ℛ∈ℂNc×M,where M=ℛ is the number of positions contained in ℛ. When there is no ambiguity about the choice of the region of interest, we shall simply use the notation Ω1.


## 2.2. Flip Angle Shimming Pulse

Letut∈ℂ1×Nc and kt∈ℝ1×3 be the RF and k-space trajectory waveforms of a given FA shimming pulse. Let T be the duration of the RF and magnetic field gradient (MFG) waveforms and ω0 (rad/s) be the resonance offset at location r due to the static field nonuniformity. Let us use below the notation Pr,t=uteikt⋅r. In the small tip angle approximation, the effective (time-averaged) control field created by this RF and magnetic field gradient (MFG) pulse pair is given by [8](3)ω¯1P,r=P¯rω1r,where(4)P¯r=1T∫t=0Tuteikt⋅reiω0rT−tdt.Note that the flip angle created by the pulseP is simply ϕP,r=Tω¯1P,r. As a result, the mean, μω¯1, and variance, σω¯12, of the effective control field ω¯1 are given by(5)μω¯1r=P¯rμr,σω¯1r2=P¯rCrP¯r∗.In this model, the effect of off-resonances is taken into account but is not modelled as a random processω0r. Equation (4) hence can be implemented as such if the static field offset is mapped in a subject-specific manner. Alternatively, one can assume ω0=0, i.e., that the static field offset has a negligible influence on pulse performance, which is acceptable for short (Tω0≪1) pulses.In the following, we present a method to compute universal FA shimming pulses which is entirely based on a second-order statistical model for the control field. Along with this statistically driven framework, we provide a method to evaluate the performance of a universal FA shimming pulse.


## 2.3. Statistical Robust Design of Universal Pulses

LetP be a candidate FA shimming pulse to create a unitary effective control field ω¯1. Let us denote by(6)δP,r=ω¯1P,r−1,the discrepancy between the unit target and ω¯1P,r, and by ΠPΔ,r the probability that δP,r exceeds a certain error level Δ>0. To optimize P, we can exploit here a statistical robust design (SRD) approach [20]. Let us introduce a new design parameter 0≤θ≤1, called SRD tolerance, and define the error level ΔP,θ as(7)ΔP,θ=minΔ≥0;ΠPΔ,r≤θ∀r∈ℛdesign,where ℛdesign denotes the region of interest for the design of the pulse (any position outside r being therefore ignored in this metric). The SRD problem may then be posed as the minimization problem:(8)minPΔP,θ.Denoting now byFδP,r the cumulative distribution function of δP,r, the condition ΠPΔ,r≤θ gives(9)FδP,rΔ≥θ.Hence, the problem described by equation (8) takes the equivalent form:(10)minPmaxr∈ℛdesignFδP,r−1θ,where FδP,r−1 denotes the inverse function of FδP,r. A useful simplification of this optimization problem is then obtained by replacing the maximum (ℒ∞ norm) by the ℒq norm for a certain choice of q≥1. This leads to the following new pulse optimization problem:(11)minP∑r∈ℛdesignFδP,r−1θq.If, now, for a given pulse candidateP, the following two conditions are satisfied for all positions r in ℛdesign:(12)varω¯1P,r≪ω¯1P,r2,ω¯1P,r−12≪varω¯1P,r,then we can show (see Appendix A) that the pulse optimization problem (11) can be further simplified to take the form:(13)minP∑r∈ℛdesignω¯1P,r−12q/2,where the dependence on the second-order B1+ statistics can be made explicit:(14)ω¯1P,r−12=12P¯rCP¯r∗+P¯rμr2−12+12P¯rCP¯r∗.Retrospectively, the link with the database-driven objective function (simultaneous design across all subjects of the database ofB1+ maps), denoted here by fUP, can be explained as follows. The latter can be defined as the mean across the design database of the mean squares error of the effective control field:(15)fUPP=∑n=1Ns∑r∈ℛdesignP¯rω1nr−12.Now, if we assimilate sample mean and statistical mean, we obtain(16)fUPP=∑r∈ℛdesignP¯rω1r−12=∑r∈ℛdesignω¯1P,r−12,which leads us to the same expression as the one obtained in equation (13) for q=2. This shows the near equivalence (note however that the objective function defined in reference [15] is proportional to the average (ℒ1 norm) of the FA-NRMSE rather than the ℒ2 norm and that the brain region for the evaluation of the FA-NRMSE was subject-specific rather than fixed as in equation (15)) of the database-driven and the statistical approaches when sample mean and statistical mean can be assimilated. This is acceptable in particular when the design database is large enough to be representative of the whole targeted population.


## 2.4. Derivation of a Lower Bound for the Coefficient of Variation of the Flip Angle

Letαr be a solution to the problem of minimizing the variance of the control field at position r among the RF shim solutions that return on average, up to a phase term, a unitary control field:(17)αr=argmina∈ℂ1×Ncaω1r−12subject toaω1r=1.In AppendixB, an explicit description of αr is obtained under the condition that(18)βr≜μrC−1rμr∗≫1.The expression forαr then reads(19)α=C−1μ∗β,and the variance of the value of the objective at the optimum is given by(20)αrCrαr∗=1β.Furthermore, it is shown in AppendixC that the locally optimal solution αr sets a performance limit for the coefficient of variation (CV) of the effective control field:(21)gr=σω¯1rμω¯1r,which reads(22)gr2≥αrCrαr∗.Hence, under the conditionβ≫1, the lower bound for the CV of the effective control field, which we may also call ultimate FA-CV in the linear STA regime, and below denoted by g∞, reads(23)g≥1β≜g∞.The conditionβ≫1, on the other hand, has a clear physical interpretation: the mean control field magnitude μ dominates its standard deviation. The validity of this condition in practice will be tested subsequently using experimental data.It is also of interest to define the ultimate FA-CV—in the sense of equation (23)—for the circularly polarized (CP) mode of the TX array. In this mode, a fixed relationship between the RF signals applied to the TX array ports is imposed in order to reproduce the behavior of a single port volume coil. The CP mode can be defined with the vector a⌣∈ℂ1×Nc which describes the way the channels of the TX array interfere. The mean and variance of the CP mode control field, μ⌣,C⌣, is then given by μ⌣=a⌣μ and C⌣=a⌣Ca⌣∗. Application of equation (23) on the CP mode control field statistics gives(24)g⌣∞=1β⌣,where(25)β⌣=a⌣μ2a⌣Ca⌣∗.In what follows, we present an application of the proposed statistical model of the intersubject variability of the control field and their implication in terms of pulse design and flip angle shimming performance assessment.


## 3. Material and Methods

### 3.1. Control Field Second-Order Statistics

The data consisted of a collection ofN=40B1+ maps measured on different adult subjects (age =40±15 years) which were recruited as part of previous studies, described in [14, 15, 21]. The MRI data were obtained on a parallel transmit-enabled whole-body Siemens 7 T system (Siemens Healthineers, Erlangen, Germany) equipped with an SC72 whole-body gradient insert (200 mT/m/ms and 100 mT/m maximum gradient slew rate and maximum gradient amplitude, respectively) and the 8Tx-32Rx Nova head coil (Nova Medical, Wilmington, MA, USA) for signal transmission and reception. This RF coil array possesses a CP mode for transmission in which the amplitude and phases at excitation ports 2 to 8 are imposed by the amplitude and phase applied to the first excitation port. This operation mode thus has only one degree of freedom (the complex amplitude applied to the first excitation port) and is designed to create the TX profile corresponding to a birdcage resonator. The parallel transmission system consisted of eight transmitters (1 kW peak power per channel) driven with the Siemens Tim Tx (Step 2.3) system. The protocol used to reconstruct the quantitative B1+ maps was acquired from a multislice interferometric turbo-FLASH acquisition [22, 23] (5 mm isotropic resolution, matrix size 40×64×40, TR = 20 s, TA = 4 min 40 s).The control field second-order statisticsμr,Cr of the TX array head coil were computed from the TX field maps of the first (chronology of data acquisition) 35 subjects of the database using the mean and covariance estimators, as described in the Theory section. For the estimation of the covariance matrix, the Ledoit–Wolf correction was used due to the limited number of subjects as compared to the dimension of the covariance matrix (64).


### 3.2. Coefficient of Variation of the Flip Angle

The lower bounds for the CV of the effective control fieldω¯1 in the linear STA regime for the CP and pTX modes were computed from the estimated second-order statistics of the control field μ^,C^ using equation (23). This enabled us in particular to compare the ultimate FA shimming performance when exploiting full pTX capability (independent RF signals applied to each port of the TX array) with the FA shimming performance of the CP mode of the TX array.


### 3.3. Statistical Robust Design ofkT-Point Pulses

For that study, the spins were assumed to be at resonanceω0=0 everywhere, i.e., we assumed that the off-resonance effect on excitation was negligible. Computed FA shimming pulses were 7-kT-point [11] pulses (7 rectangular subpulses), wherein the duration of the RF subpulses and MFG blips were set to 140 μs and 50 μs, respectively. In all optimizations, the target flip angle, ϕtarget, was 10 deg. Furthermore, the RF peak power applied at each excitation port was limited to 500 W, the total RF energy (respectively, the RF energy per TX channel) of the pulse was assigned to not exceed 60 mJ (respectively, 12 mJ), and the k-space displacement in x, y, and z created by each MFG blip was imposed to not exceed 50 rad/m.A first FA shimming pulse, referred to as SRD-UP, was designed from the second-order control field statistics of the Nova TX array head coil (estimated inSection 3.1) by solving numerically the problem (same as equation (13) but nonunitary target, T duration of the pulse):(26)minP∑r∈ℛdesignω¯1P,r−ϕtargetT2q/2.In the above expression,ℛdesign was defined as the set of positions that belonged to the brain region in at least half of the subjects of the database, i.e.,(27)r∈ℛdesign⟺1N∑n=1N1ℛbrainnr≥0.5,where ℛbrainn denotes the mask of the brain for the n-th subject of the database and 1ℛbrainn⋅ is the indicator function of this region.For comparison with SRD-UP, the solution (same parameterization and constraints as for the design of SRD-UP) of the database-driven UP objective was also computed. The database-driven UP objective was defined as follows:(28)minP∑n=1Nsℛbrainn−1∑r∈ℛbrainnP¯rω1nr−ϕtargetT2.This pulse is referred to as UP. Finally, a standard hard pulse where the transmit array is used in CP mode was defined. The amplitude of this pulse thus was set to return on average on the design database the target flip angle.We then studied the spatial distribution of the mean and CV of the flip angle for UP, SRD-UP, and CP and compared it with the ultimate FA shimming performanceg∞. For convenience, this performance is attributed to a virtual pulse P which would satisfy (see equations (4) and (19)) TP¯r=ϕtargetαr. In the following, this pulse is referred to as UP∞. Note that αr being a local solution, no RF and MFG waveforms can satisfy strictly this relationship unless the k-space trajectory defined by this pulse offers a complete coverage of the transmit k-space trajectory. This would lead in practice to a prohibitively long FA shimming pulse.Using the five last subjects of the cohort (the 35 first constituting the design database), the flip angle normalized root mean squares error (FA-NRMSE) of UP, SRD-UP, and CP (P one of these pulses, ϕBlochP,r, the FA created by pulse P at location r, 36≤n≤40):(29)1ℛbrainn∑r∈ℛbrainnϕBlochP,rϕtarget−121/2,was computed on each subject. For the FA-NRMSE computation, the FA maps, ϕBlochP,⋅, were obtained by numerical integration of Bloch’s equations (whereby T1 and T2 relaxations were neglected) rather than using the linear small tip angle approximation (equation (3) For every subject of this test database, a subject-tailored kT-point pulse (minimization of the nth term of the sum in equation (28)), possessing the same parameterization as the above universal pulses (SRD-UP and UP), was also designed to minimize the FA-NRMSE. These solutions are referred to as ST36, …, ST40.All pulse optimizations were performed using the second-order active-set algorithm, available in the Optimization Toolbox of MATLAB R2019a (The Mathworks, Natick, MA), which allows solving constrained optimization problems numerically. The active-set algorithm used the analytical expressions for the Jacobian of the objective function and the RF power, RF energy, and MFG constraints (Hoyos-Idrobo et al. [24]). The number of iterations was limited to 2×103. The algorithm however could stop before exceeding this number of iterations provided it reached a local minimum within default tolerance.


## 3.1. Control Field Second-Order Statistics

The data consisted of a collection ofN=40B1+ maps measured on different adult subjects (age =40±15 years) which were recruited as part of previous studies, described in [14, 15, 21]. The MRI data were obtained on a parallel transmit-enabled whole-body Siemens 7 T system (Siemens Healthineers, Erlangen, Germany) equipped with an SC72 whole-body gradient insert (200 mT/m/ms and 100 mT/m maximum gradient slew rate and maximum gradient amplitude, respectively) and the 8Tx-32Rx Nova head coil (Nova Medical, Wilmington, MA, USA) for signal transmission and reception. This RF coil array possesses a CP mode for transmission in which the amplitude and phases at excitation ports 2 to 8 are imposed by the amplitude and phase applied to the first excitation port. This operation mode thus has only one degree of freedom (the complex amplitude applied to the first excitation port) and is designed to create the TX profile corresponding to a birdcage resonator. The parallel transmission system consisted of eight transmitters (1 kW peak power per channel) driven with the Siemens Tim Tx (Step 2.3) system. The protocol used to reconstruct the quantitative B1+ maps was acquired from a multislice interferometric turbo-FLASH acquisition [22, 23] (5 mm isotropic resolution, matrix size 40×64×40, TR = 20 s, TA = 4 min 40 s).The control field second-order statisticsμr,Cr of the TX array head coil were computed from the TX field maps of the first (chronology of data acquisition) 35 subjects of the database using the mean and covariance estimators, as described in the Theory section. For the estimation of the covariance matrix, the Ledoit–Wolf correction was used due to the limited number of subjects as compared to the dimension of the covariance matrix (64).


## 3.2. Coefficient of Variation of the Flip Angle

The lower bounds for the CV of the effective control fieldω¯1 in the linear STA regime for the CP and pTX modes were computed from the estimated second-order statistics of the control field μ^,C^ using equation (23). This enabled us in particular to compare the ultimate FA shimming performance when exploiting full pTX capability (independent RF signals applied to each port of the TX array) with the FA shimming performance of the CP mode of the TX array.


## 3.3. Statistical Robust Design ofkT-Point Pulses

For that study, the spins were assumed to be at resonanceω0=0 everywhere, i.e., we assumed that the off-resonance effect on excitation was negligible. Computed FA shimming pulses were 7-kT-point [11] pulses (7 rectangular subpulses), wherein the duration of the RF subpulses and MFG blips were set to 140 μs and 50 μs, respectively. In all optimizations, the target flip angle, ϕtarget, was 10 deg. Furthermore, the RF peak power applied at each excitation port was limited to 500 W, the total RF energy (respectively, the RF energy per TX channel) of the pulse was assigned to not exceed 60 mJ (respectively, 12 mJ), and the k-space displacement in x, y, and z created by each MFG blip was imposed to not exceed 50 rad/m.A first FA shimming pulse, referred to as SRD-UP, was designed from the second-order control field statistics of the Nova TX array head coil (estimated inSection 3.1) by solving numerically the problem (same as equation (13) but nonunitary target, T duration of the pulse):(26)minP∑r∈ℛdesignω¯1P,r−ϕtargetT2q/2.In the above expression,ℛdesign was defined as the set of positions that belonged to the brain region in at least half of the subjects of the database, i.e.,(27)r∈ℛdesign⟺1N∑n=1N1ℛbrainnr≥0.5,where ℛbrainn denotes the mask of the brain for the n-th subject of the database and 1ℛbrainn⋅ is the indicator function of this region.For comparison with SRD-UP, the solution (same parameterization and constraints as for the design of SRD-UP) of the database-driven UP objective was also computed. The database-driven UP objective was defined as follows:(28)minP∑n=1Nsℛbrainn−1∑r∈ℛbrainnP¯rω1nr−ϕtargetT2.This pulse is referred to as UP. Finally, a standard hard pulse where the transmit array is used in CP mode was defined. The amplitude of this pulse thus was set to return on average on the design database the target flip angle.We then studied the spatial distribution of the mean and CV of the flip angle for UP, SRD-UP, and CP and compared it with the ultimate FA shimming performanceg∞. For convenience, this performance is attributed to a virtual pulse P which would satisfy (see equations (4) and (19)) TP¯r=ϕtargetαr. In the following, this pulse is referred to as UP∞. Note that αr being a local solution, no RF and MFG waveforms can satisfy strictly this relationship unless the k-space trajectory defined by this pulse offers a complete coverage of the transmit k-space trajectory. This would lead in practice to a prohibitively long FA shimming pulse.Using the five last subjects of the cohort (the 35 first constituting the design database), the flip angle normalized root mean squares error (FA-NRMSE) of UP, SRD-UP, and CP (P one of these pulses, ϕBlochP,r, the FA created by pulse P at location r, 36≤n≤40):(29)1ℛbrainn∑r∈ℛbrainnϕBlochP,rϕtarget−121/2,was computed on each subject. For the FA-NRMSE computation, the FA maps, ϕBlochP,⋅, were obtained by numerical integration of Bloch’s equations (whereby T1 and T2 relaxations were neglected) rather than using the linear small tip angle approximation (equation (3) For every subject of this test database, a subject-tailored kT-point pulse (minimization of the nth term of the sum in equation (28)), possessing the same parameterization as the above universal pulses (SRD-UP and UP), was also designed to minimize the FA-NRMSE. These solutions are referred to as ST36, …, ST40.All pulse optimizations were performed using the second-order active-set algorithm, available in the Optimization Toolbox of MATLAB R2019a (The Mathworks, Natick, MA), which allows solving constrained optimization problems numerically. The active-set algorithm used the analytical expressions for the Jacobian of the objective function and the RF power, RF energy, and MFG constraints (Hoyos-Idrobo et al. [24]). The number of iterations was limited to 2×103. The algorithm however could stop before exceeding this number of iterations provided it reached a local minimum within default tolerance.


## 4. Results

### 4.1. Second-Order Field Statistics

In Figure1, the second-order statistics of the control field for the Nova TX array are displayed for one axial slice through the brain. The mean control field μ thus gives rise to eight axial images, each one corresponding to one TX channel. We also provide the mean control field of the CP mode, μ⌣=a⌣μ. The covariance matrix is composed of the control field variance of every TX channel (eight diagonal images) as well as the covariance between channels (off-diagonal images).Figure 1
(a) Axial view of the mean control field (in Hz/volt) of the Nova TX array (images 1 to 8 correspond to the 8 TX channels and the last image correspond to the CP excitation mode) and (b) the covariance of the control field (Hz2/volt2) for the Nova TX array.
(a)(b)


### 4.2. Ultimate FA Shimming Performance

Figure2 displays sagittal, coronal, and axial images of the lower bound for the CV of the effective control field (ultimate FA shimming performance) in CP mode and the pTX operation mode. Here, it can be seen that g∞<g⌣∞, i.e., that the pTX drive mode allows for a significant decrease in the CV of the effective control field as compared to the CP mode. In other words, for every voxel in the brain, there exists a pTX excitation mode that yields less intersubject variability than the CP mode. Quantitative comparisons of these performance boundaries are also provided in the form of cumulative histograms of g∞ and g⌣∞ in Figure 3. Finally, histograms of the distributions of the ratio of CVs g⌣∞/g∞ are provided in Figure 4. From this plot, it can be seen in particular that the CV of the effective control field ω¯1 can be reduced by a factor greater than two in some locations in the brain with pTX as compared to the CP mode.Figure 2
Sagittal, coronal, and axial views of the lower bound for CV of the effective control field (ω¯1) for the CP drive mode (left), g⌣∞, and the full pTX drive mode (right), g∞. As compared to the CP drive mode, FA shimming performances, measured as the CV of the effective control field on each location, are clearly superior using pTX (i.e., 8 degrees of freedom as opposed to 1 for the CP mode). It also shows that in CP mode, the greatest intersubject variability occurs in the cerebellum and in the temporal lobes.Figure 3
Cumulative histograms of the lower bound for the control field CV for the CP (dashed curves) and pTX (solid curves) drive modes. With pTX, the 80% quantile of control field CV is<9%. Without pTX (CP mode), it increases to ∼15%.Figure 4
Cumulative histogram of the ratiog⌣∞/g∞, whereby g∞ and g⌣∞ are the theoretical lower bounds for the effective control’s CV in the pTX (8 degrees of freedom) and the CP (1 degree of freedom) modes. The more “agile” pTX drive mode can adapt to the known control field statistics and thereby allows reducing significantly the lower bound for the CV of the effective control field. A reduction factor of at least 1.5 is obtained for 80% of the voxels of the brain. In some areas, the CV reduction factor reaches a factor of 2 or more.


### 4.3. FA Shimming Performance Simulation

The FA shimming performance analysis ofUP∞, UP, SRD-UP, and for the standard hard pulse in the CP mode is shown in Figures 5 and 6. For UP and SRD-UP, in the brain (in the sense of (30)), the mean flip angle ϕr (statistical mean) is very close to the target (10deg) as it remains in the interval 9,11deg. By definition of α (see equation (19)), for the virtual pulse UP∞, the mean flip angle is exactly equal to the target. As it can be expected, for the CP mode, the mean flip angle can deviate very strongly from the 10deg target as it reaches for instance in the center >20deg. On the other hand, the median CV of the flip angle amounts to 12% for UP and CP and 10% for SRD-UP and is ultimately <8% for UP∞. In the “worst location” which was found to be the cerebellum (see Figure 5), the CV of the flip angle does not exceed 15% with UP∞ but can exceed 20% with “real” pulses.Figure 5
Median sagittal view of the mean (a) and (b) CV of the flip angle for (from top to bottom)UP∞, UP, SRD-UP, and CP. As expected, for all pulses except CP, the mean flip angle FA is very close to the 10deg target. For the CV of the flip angle, it is of interest to note the gain in the performance of SRD-UP as compared to CP (consistent with the result of Figure 2).
(a)(b)Figure 6
Mean value (a, b) and (c) CV of the flip angle forUP∞, UP, SRD-UP, and CP. For each pulse, the minimum, the maximum, the 75%, and the 50% quantiles of the flip angle mean and CV values are represented. For better readability, due to the large variability in the mean flip angle for the CP mode, the corresponding bar plot was displayed in a separate graph (b). Note that UP∞ returns by definition an unbiased flip angle (hence the absence of box for this nonrealizable pulse).
(a)(b)(c)The result of the FA-NRMSE analysis made separately on the design and test subjects is presented in Figure7. The results obtained on the test subjects are compliant with the FA-NRMSE statistics obtained on the design database in the sense that no significant offset (unpaired Wilcoxon rank-sum test) was found between the FA-NRMSE obtained in the design and test databases. On the test subjects, the performance was very comparable between UP∞, UP, and SRD-UP, with a median FA-NRMSE of 5 to 7%. The subject-tailored pulses returned FA-NRMSEs <2%.Figure 7
FA-NRMSE performances of (a)UP∞, UP, SRD-UP, and ST36, …, ST40 (ST) and (b) CP on the design database (box plots indicating the minimum, maximum, 25%, 50% (red horizontal line), and 75% quantiles of the FA-NRMSE values) and the test database (5 subjects) (blue stars). Note the absence of box plot for ST as these pulses are defined only for the five test subjects. For the FA-NRMSE performance assessment on the test database, individual performances are represented with blue stars and median FA-NRMSE across the five subjects is represented with the dashed blue horizontal line.
(a)(b)


### 4.4. Storage and Computation Time Efficiency

Denoting byNv the average number of voxels in the brain, the storage volume for the first- and second-order statistics of the control field and for a database of N subjects are Nc×Nc+3×Nv/2 and N×Nv×Nc, respectively. Thus, for N≥Nc+3/2 (6 in our case), the statistical representation becomes more parsimonious. From the computation time point of view, the processing time (DELL Precision 7540, 32 Gb of RAM, 8 cores, Intel CORE i9 vPro 9th Gen) was 3270 s (669 iteration and 1851 evaluations of the objective function) and 160 s (489 iterations and 1471 evaluations of the objective function) for the computation of UP and SRD-UP, respectively. The computation time for one evaluation of the objective function (and the Jacobian) was 2.13 s and 0.26 s for UP and SRD-UP, respectively. As these workloads scaled linearly with the number of voxels (4.22×105 in total for the 35 subjects of the database and 1.1×104 for ℛdesign (see equation (27))), we found that the workloads per voxel were 5 μs for UP and 17.4 μs for the computation of UP and SRD-UP, respectively. As a result, but keeping in mind that these numbers may depend on the pulse parameterization (e.g., number of kT-points), the design of SRD-UP is computationally more efficient than the UP design when the size of the design database is greater than 17.4/5≃3.5.


## 4.1. Second-Order Field Statistics

In Figure1, the second-order statistics of the control field for the Nova TX array are displayed for one axial slice through the brain. The mean control field μ thus gives rise to eight axial images, each one corresponding to one TX channel. We also provide the mean control field of the CP mode, μ⌣=a⌣μ. The covariance matrix is composed of the control field variance of every TX channel (eight diagonal images) as well as the covariance between channels (off-diagonal images).Figure 1
(a) Axial view of the mean control field (in Hz/volt) of the Nova TX array (images 1 to 8 correspond to the 8 TX channels and the last image correspond to the CP excitation mode) and (b) the covariance of the control field (Hz2/volt2) for the Nova TX array.
(a)(b)


## 4.2. Ultimate FA Shimming Performance

Figure2 displays sagittal, coronal, and axial images of the lower bound for the CV of the effective control field (ultimate FA shimming performance) in CP mode and the pTX operation mode. Here, it can be seen that g∞<g⌣∞, i.e., that the pTX drive mode allows for a significant decrease in the CV of the effective control field as compared to the CP mode. In other words, for every voxel in the brain, there exists a pTX excitation mode that yields less intersubject variability than the CP mode. Quantitative comparisons of these performance boundaries are also provided in the form of cumulative histograms of g∞ and g⌣∞ in Figure 3. Finally, histograms of the distributions of the ratio of CVs g⌣∞/g∞ are provided in Figure 4. From this plot, it can be seen in particular that the CV of the effective control field ω¯1 can be reduced by a factor greater than two in some locations in the brain with pTX as compared to the CP mode.Figure 2
Sagittal, coronal, and axial views of the lower bound for CV of the effective control field (ω¯1) for the CP drive mode (left), g⌣∞, and the full pTX drive mode (right), g∞. As compared to the CP drive mode, FA shimming performances, measured as the CV of the effective control field on each location, are clearly superior using pTX (i.e., 8 degrees of freedom as opposed to 1 for the CP mode). It also shows that in CP mode, the greatest intersubject variability occurs in the cerebellum and in the temporal lobes.Figure 3
Cumulative histograms of the lower bound for the control field CV for the CP (dashed curves) and pTX (solid curves) drive modes. With pTX, the 80% quantile of control field CV is<9%. Without pTX (CP mode), it increases to ∼15%.Figure 4
Cumulative histogram of the ratiog⌣∞/g∞, whereby g∞ and g⌣∞ are the theoretical lower bounds for the effective control’s CV in the pTX (8 degrees of freedom) and the CP (1 degree of freedom) modes. The more “agile” pTX drive mode can adapt to the known control field statistics and thereby allows reducing significantly the lower bound for the CV of the effective control field. A reduction factor of at least 1.5 is obtained for 80% of the voxels of the brain. In some areas, the CV reduction factor reaches a factor of 2 or more.


## 4.3. FA Shimming Performance Simulation

The FA shimming performance analysis ofUP∞, UP, SRD-UP, and for the standard hard pulse in the CP mode is shown in Figures 5 and 6. For UP and SRD-UP, in the brain (in the sense of (30)), the mean flip angle ϕr (statistical mean) is very close to the target (10deg) as it remains in the interval 9,11deg. By definition of α (see equation (19)), for the virtual pulse UP∞, the mean flip angle is exactly equal to the target. As it can be expected, for the CP mode, the mean flip angle can deviate very strongly from the 10deg target as it reaches for instance in the center >20deg. On the other hand, the median CV of the flip angle amounts to 12% for UP and CP and 10% for SRD-UP and is ultimately <8% for UP∞. In the “worst location” which was found to be the cerebellum (see Figure 5), the CV of the flip angle does not exceed 15% with UP∞ but can exceed 20% with “real” pulses.Figure 5
Median sagittal view of the mean (a) and (b) CV of the flip angle for (from top to bottom)UP∞, UP, SRD-UP, and CP. As expected, for all pulses except CP, the mean flip angle FA is very close to the 10deg target. For the CV of the flip angle, it is of interest to note the gain in the performance of SRD-UP as compared to CP (consistent with the result of Figure 2).
(a)(b)Figure 6
Mean value (a, b) and (c) CV of the flip angle forUP∞, UP, SRD-UP, and CP. For each pulse, the minimum, the maximum, the 75%, and the 50% quantiles of the flip angle mean and CV values are represented. For better readability, due to the large variability in the mean flip angle for the CP mode, the corresponding bar plot was displayed in a separate graph (b). Note that UP∞ returns by definition an unbiased flip angle (hence the absence of box for this nonrealizable pulse).
(a)(b)(c)The result of the FA-NRMSE analysis made separately on the design and test subjects is presented in Figure7. The results obtained on the test subjects are compliant with the FA-NRMSE statistics obtained on the design database in the sense that no significant offset (unpaired Wilcoxon rank-sum test) was found between the FA-NRMSE obtained in the design and test databases. On the test subjects, the performance was very comparable between UP∞, UP, and SRD-UP, with a median FA-NRMSE of 5 to 7%. The subject-tailored pulses returned FA-NRMSEs <2%.Figure 7
FA-NRMSE performances of (a)UP∞, UP, SRD-UP, and ST36, …, ST40 (ST) and (b) CP on the design database (box plots indicating the minimum, maximum, 25%, 50% (red horizontal line), and 75% quantiles of the FA-NRMSE values) and the test database (5 subjects) (blue stars). Note the absence of box plot for ST as these pulses are defined only for the five test subjects. For the FA-NRMSE performance assessment on the test database, individual performances are represented with blue stars and median FA-NRMSE across the five subjects is represented with the dashed blue horizontal line.
(a)(b)


## 4.4. Storage and Computation Time Efficiency

Denoting byNv the average number of voxels in the brain, the storage volume for the first- and second-order statistics of the control field and for a database of N subjects are Nc×Nc+3×Nv/2 and N×Nv×Nc, respectively. Thus, for N≥Nc+3/2 (6 in our case), the statistical representation becomes more parsimonious. From the computation time point of view, the processing time (DELL Precision 7540, 32 Gb of RAM, 8 cores, Intel CORE i9 vPro 9th Gen) was 3270 s (669 iteration and 1851 evaluations of the objective function) and 160 s (489 iterations and 1471 evaluations of the objective function) for the computation of UP and SRD-UP, respectively. The computation time for one evaluation of the objective function (and the Jacobian) was 2.13 s and 0.26 s for UP and SRD-UP, respectively. As these workloads scaled linearly with the number of voxels (4.22×105 in total for the 35 subjects of the database and 1.1×104 for ℛdesign (see equation (27))), we found that the workloads per voxel were 5 μs for UP and 17.4 μs for the computation of UP and SRD-UP, respectively. As a result, but keeping in mind that these numbers may depend on the pulse parameterization (e.g., number of kT-points), the design of SRD-UP is computationally more efficient than the UP design when the size of the design database is greater than 17.4/5≃3.5.


## 5. Discussion

In this work, we have shown that the statistical robust and the UP design methods are theoretically nearly equivalent when the size of the database is sufficiently large to embrace the variability of the resonant TX RF field across the population of interest. Despite its apparent complexity, the statistical robust design algorithm defined by equation (13) has a rather straightforward implementation as the general term of the summation across voxels depends explicitly on the mean vector and covariance matrix of the TX RF field (see equation (14)). As compared to the database-driven UP design methodology (Gras et al. [14]), an important difference is the involvement of the covariance matrix Cr which has not equivalent in the original formulation. To our knowledge, the design of UPs using the SRD formalism, originally applied by Alotto and colleagues for the design of electromagnetic structures, has not been proposed so far for RF pulse design and could represent some advantages. First, as the SRD computation time does not depend on the size (number or subjects) of the database, this approach enables integration of potentially very large databases, which is a priori necessary to ensure optimal UP performance. Secondly, equipped with this model, one can likewise compute probabilities for obtaining a certain deviation from a target flip angle. Doing such estimations using a database of TX field maps would require unrealistically large number of subjects.The theory presented in this work defines the performance of a FA shimming pulse as the coefficient of variation (std/mean) of the flip angle. Furthermore, in the linear STA regime, it is shown that this CV cannot be lower than a certain valueg∞ which depends on the local second-order statistics of the control field. In the process of validating—and ultimately certifying—the use of pTX in a clinical setup, this performance measure may prove useful. One concrete application of this concept is the comparison of the CP mode ultimate performance (g⌣∞) with the pTX performance (g∞): for the Nova TX array, we found that the availability of eight independent TX channels to locally tune the drive mode so as to minimize the variance of the effective control allowed reducing the CV of the flip angle by a factor of ≃2, while achieving close to the target on average (see Figure 4). Interestingly, this theoretical analysis sheds some light onto the concept of UPs. Indeed, with statistical arguments, one can take advantage of the TX array architecture to mitigate both the spatial heterogeneity and the intersubject variability of the transmit RF field in the brain. Finally, the coefficient of variation metric for UP∞ (see equation (23)), being characteristic of the TX array architecture, can potentially serve as a guide for coil-design purposes.The present study was performed with theB1+ data collected on 40 subjects, from which a subgroup of 35 subjects was used to compute statistics. Given the high dimension of the covariance matrix of the control field for an eightfold TX array system, it was necessary to use the Ledoit–Wolf correction for the estimation of the covariance matrix of the control field (to obtain positive definite estimates). However, we can still expect a bias in the estimation of the covariance matrix. As a result, it is possible that the CV analysis reported here (lower bound and CV maps for the different kT-point pulses computed in this study) can deviate from the analysis that one would obtain with a larger database to estimate the control field second-order statistics.The theory developed in this study is based on the assumption that the statistical distribution of the control field vectorω1r is Gaussian. This assumption is very convenient from the mathematical point of view as this property is inherited by the effective control field ω¯1r in the linear STA regime. With a few tens of subjects, this hypothesis has not yet been verified quantitatively, but given the complexity and the multiplicity of effects leading to the observed variability, the Gaussian model is a reasonable assumption a priori. Another hypothesis made in this study is the fact that β=μ∗C−1μ is much greater than one. Here we notice that this is equivalent to the assumption that g∞ is much smaller than one. In Figure 2, g∞ was seen not to exceed 0.2, confirming that this hypothesis was approximately valid over the whole brain.Extending this work to large flip angle pulse is one interesting prospect. However, in this case, the linear relationship betweenω1 and ω¯1 exploited to derive a simple expression of the SRD objective does not hold. Here a pragmatic approach would be to replace exact calculation of the first and second moments of the flip angle (seen as a random variable) by numerical approximations thereof. A possible direction to perform such calculation is the use of the so-called unscented transform method [25] which seems here particularly adapted [20].An important difference with the pulse design method reported in [13] is the deterministic nature of the off-resonance term ω0 in the current statistical design. In this study, we set ω0=0 for the design of UPs using the SRD formalism, i.e., we neglected off-resonance effects in the performance assessment. While acceptable for short excitation pulses, this can be a significant source of error in practice when the pulse duration becomes large [26]. Extending the statistical robust design approach to account for a random off-resonance term seems however feasible with some additional hypotheses. Hence, if it is assumed that ω0r is a centered Gaussian random variable independent of ω1, it remains possible in equation (3) to compute the mean deviation of the effective control field with respect to the effective control field at resonance. This quantity takes the form:(30)ω¯1r−ω¯10r21/2,where ω¯10 denotes the effective control field at resonance, and where the expectation operator ⋅ covers ω1 as well as ω0 statistics. A possible direction is to exploit this quantity as an additional objective to improve the robustness against off-resonances while taking into account the regional dependence of the variance of the resonance frequency offset.In this work, for simplicity, the statistical robust pulse design implementation did not include specific absorption rate (SAR) constraints. However, as it was done in the past for the design of UPs, there would be no obstacle in adding explicit local SAR constraint to the statistical robust design procedure in the form of virtual observation points [27]. Interestingly also, the RF electric field distribution (E-field) at the origin of the SAR is another imperfectly controlled physical quantity affecting RF safety management procedure, in particular at UHF and with transmit array architectures [28]. Like the B1+ field, the E-field could advantageously be treated as well as a stochastic process, thereby allowing for a probabilistic model of the local 10 g SAR. Naturally, the dual problem consisting of minimizing the SAR under performance (flip angle error) constraints [29] could benefit as well from such a unified statistical model. One obstacle for the estimation of the E-field statistics yet remains the absence of direct method to measure the E-field in vivo and thus the necessity to resort to electromagnetic simulations on a large amount of head models to give access to them. Furthermore, due to the interweaved nature of the RF electric and magnetic fields, the values of the E- and B1+-field, seen here as random processes, are in fact necessarily correlated.


## 6. Conclusions

We have presented in this work a statistical robust approach to design calibration-free pTX universal pulses for brain imaging at ultrahigh field. Experimentally, this method has been shown to provide comparable FA-NRMSE results as the conventional database-driven UP design, but presents a significant advantage in terms of computational efficiency when the size of the design database exceeds 4. Furthermore, in contrast with the database-driven design of UP, the computation time of the statistical robust UP does not depend on the size of this database and consequently enables integration of a very large design database. This is an important criterion to ensure optimal UP performance. Furthermore, the theoretical lower bound for the coefficient of variation of the flip angle derived in this work,g∞, appears as a valuable analytical performance limit to assess the FA shimming performance of a given pulse. Finally, this metric can be useful to characterize, and possibly also optimize, the performance of TX array architectures.



---

*Source: 6018107-2020-09-18.xml*