Design of a shim coil array matched to the human brain anatomy

The purpose of this study is to introduce a novel design method of a shim coil array specifically optimized for whole brain shimming and to compare the performance of the resulting coils to conventional spherical harmonic shimming.

MR frequency variations scale linearly with the main magnetic field.
In order to improve the field homogeneity, conventional second-order spherical harmonic (SH) shim coils are used in routine scanners. Typically, prior to execution of imaging or spectroscopy sequences, the input currents of the shim coils are adjusted to achieve a certain reduction of the field inhomogeneity. However, even after the best possible secondorder SH shimming, residual B 0 inhomogeneity still exists, especially in local regions in the vicinity of air cavities and canals such as frontal lobes or temporal lobes of the human brain. Currently, high-order spherical harmonics coils, 5,6 multi-coils, [7][8][9][10][11] and integrated RF/B 0 shim coil arrays [12][13][14] have been proposed and demonstrated to provide better shimming capabilities than the standard manufacturer-provided shimming. The mentioned high-order SH coils were designed to generate magnetic fields with an orthonormal set of pure SH shapes. By contrast, the elements of multi-coil shim arrays and integrated RF/B 0 shim arrays typically have a simple fixed shape, such as a circular loop and are distributed in a regular pattern over a surface such as a cylinder or closefitting helmet. During the design of all these generic coil arrays and their elements, typical B 0 field maps acquired in human subjects have not been considered explicitly.
In order to reduce the field inhomogeneity over a specific region of the human brain, such as the medial temporal lobe, a single coil shim insert 15 has been proposed and designed for magnetic resonance spectroscopy. During the design, the average of field maps of 10 normal human subjects has been used as a target field over a spherical region. Such an approach, however, is not applicable to whole brain shimming as a fixed mask cannot be applied to each of the 10 subjects since individual brains have different geometries. Moreover, the usage of a fixed mask tends to exclude regions with high inhomogeneities, typically located close to the periphery.
If custom shim coils with more than one channel are designed, averaging is not sufficient. Therefore, alternative approaches 16,17 were proposed that consider field maps of multiple human subjects simultaneously to find a basis set that represents the particular spatial patterns of field inhomogeneity occurring in the brain using substantially fewer basis field terms than the number of input field maps. In particular, the approach in Li et al 17 relies on a singular value decomposition (SVD) applied to multiple field maps of the brain and uses the singular vectors as shim basis fields. To account for individual shape differences, brain normalization 18 was applied prior to the SVD. The required basis fields could provide better B 0 homogeneity than the first through third-order SH fields. However, the derived basis field patterns may be challenging to generate by a shim coil in practice and may not even satisfy Maxwell's equations. Furthermore, the validity of the spatial normalization step is questionable. On the one hand, it applies spatial transformations not compatible with Maxwell's equations. On the other hand, it removes rigid head motion, which makes one of the major factors affecting the susceptibility-induced frequency distribution invisible for the coil design step.
In this work, a design method of a shim coil array specifically optimized for whole-brain shimming in humans is proposed. The shim design is based on the stream function method 19 applied to 84 target field maps acquired in 12 normal volunteers in seven different head positions, followed by the singular value decomposition. The leave-one-out cross validation technique, commonly used in machine learning, 20 is applied here to determine the total number of required elements of the coil array. Statistical inference based on a bootstrap approach 21 is introduced to verify the stability of the coil array design. Numerical results are presented to demonstrate the validity of the proposed design approach.

| METHODS
The proposed design method is summarized as follows. First, the field maps of multiple human brains are acquired and preprocessed. The resulting field maps are specified as target fields and a stream function of the electric current density 19 is designed for every field map on a pre-defined currentcarrying surface. In this way, for every human brain in every position, a separate optimal shim coil is calculated. All coils share the same current-carrying surface. As the currentcarrying surface and its mesh remain the same during the design of all coil layouts, this initial design step maps each individual field map with its specific spatial position and the masked region into the common space of the electric current density on the current-carrying surface. In the following step, the matrix consisting of all resulting stream functions is analyzed by an SVD, and the left singular vectors with the largest singular values are used to obtain the wire patterns of elements of the shim coil array. A leave-one-out cross validation technique is applied to find a suitable number of channels and to verify the generalization ability of the coil array. The bootstrap approach is then used to assess whether the current sample of brain field maps is sufficient and estimate the appropriate sample size for future experiments.

| Measurement and preprocessing of human brain field maps
In this study, anatomical images and B 0 field maps of heads of 18 different healthy adult volunteers (age: 30.8 ± 3.97, gender: 14 males and 4 females) were measured on a 3T scanner (Siemens Healthcare, Erlangen, Germany) by a field mapping sequence based on a multi-echo 3D gradient echo (GRE) sequence. Before these measurements, the heads in the central position were shimmed to the second order using the routines provided by the scanner. The imaging parameters were TE = 3.84/6.7/10.56 ms, TR = 14 ms, flip angle = 8 • , field of view (FOV) = (256 × 256 × 224) mm3, and 128 × 128 × 112 matrix size. The measured maps have an isotropic resolution of 2 mm. The selected FOV was large enough to cover the whole brain (including: cerebellum, midbrain, and pons) of all subjects. All field maps from the 18 subjects were divided into two sets: the field maps from 12 subjects were designated as a training set for designing the shim coil array, while the field maps from the remaining six subjects constituted a test set. The field maps in the test set were only used to assess the performance of the final designed coil array in the following Section 3.4. In other subsections, only the field maps from the training set have been used.
In order to assess the impact of head motion on B 0 maps, each subject was measured seven times in different head positions. The first field map was measured in a reference head position considered optimal by an experienced radiographer. The other six maps were measured after rotation along x, y, and z axes by approximately ±5 • . Here, a 2D multislice echoplanar imaging (EPI) sequence was used to quantify the head movement using a prospective acquisition correction algorithm, 22 and the subjects were guided via the intercom system to comply with the requirements of the protocol.
The non-brain areas of all field maps were removed using the BET (Brain Extraction Tool) 23 based on the magnitude images from the scanned field maps at the second echo. The field maps were then unwrapped using the phase unwrapping algorithm in Robinson et al 24 based on multi-echo phase images. Magnitude of the three-dimensional gradient of the phase of the unwrapped maps weighted by the magnitude images was used to detect outliers (such as noise pixels at the boundary and occasional residual phase wraps) and modify the usable brain mask. Disconnected regions were removed from the updated mask. No further topologic operations (such as erosion) were applied to the mask to maximally preserve reliably detected field inhomogeneities close to the surface of the brain. A weak 7-voxel median filter constrained to the mask was applied to the unwrapped field maps to remove the residual outliers. Histograms of selected field maps (e.g. with the highest standard deviations within the data set) were visually inspected to confirm the reliable mask generation and the absence of apparent outliers. The processed field maps were downsampled to a 56 × 64 × 64 matrix in order to reduce the computational cost during the shim coil design step. Considering that in modern scanners the subject can be moved along the z-axis of the magnet to the optimal position without manual subject repositioning by adjusting the bed position, the median of the z-coordinates of points in the brain mask was translated to the magnet iso-center along the z-axis. Moreover, since clinical MRI scanners can rapidly adjust a radiofrequency (RF) synthesizer to compensate zeroth-order B 0 offsets 25 and gradient offsets for first-order shimming without any relevant constraints, the translated and masked brain field maps were preprocessed to remove any zeroth-and first-order SH contributions. The residual field inhomogeneities in the whole brain including midbrain, pons, and cerebellum were used for all design and validation tasks. The mean of the rootmean-square (RMS) of the residual field calculated across all 126 data sets was 26.1 ± 3.27 Hz (mean ± standard deviation).

| Design method of a shim coil for a specific brain field map
In order to obtain coil arrays for shimming all field maps in the training set, a wire pattern for each individual field map was first designed using the stream function method. The field map with the reversed polarity was used as a target magnetic field. A stream function on a cylindrical surface of certain dimensions (see below) was then optimized to minimize the standard deviation (SD) of the residual magnetic field over the whole brain. During the optimization procedure, the power dissipated by the coil was constrained in order to obtain a physically and technically meaningful solution. Based on the above considerations, the following optimization problem was used to design a shim coil for a field map: Here, ψ denotes the scalar piecewise-linear stream function 19 of the surface electric current density vector ⃗ J( ) where on a current-carrying surface Γ with a normal unit vector ⃗ n. According to Peeren,19 ψ is defined as the line integral of the cross product of ⃗ n and ⃗ J along a curve on Γ and (⃗ x) indicates the current through the curve from a fixed point to ⃗ x. B z is the z-component of the magnetic field ⃗ B generated by current densities ⃗ J on Γ, which is calculated using the Biot-Savart law. B m z indicates the measured in vivo ΔB 0 field map. The points ⃗ x i = (x i , y i , z i ) T , i = 1, …, m denote the coordinate vectors of m test points in the ROI, τ indicates the thickness of the surface Γ, σ is the electrical conductivity of the surface Γ, and P T max denotes the maximum power dissipated by the coil.
In order to house a compact transmit-receive RF head coil assembly used in previous studies, 26 a cylindrical surface with a radius of 180 mm and a length of 300 mm was used to design shim coils in this work. The current-carrying surface has a thickness of 1 mm and its material property of copper has an electrical conductivity of 5.998 × 10 7 S/m. All optimization problems were solved with two routines cgsvd and lsqi in the MATLAB (The MathWorks, Natick, Massachusetts) package named REGULARIZATION TOOLS. 27 (1) The solution of the optimization problem (1) depends on the selection of P T max . Selecting exceedingly high P T max values may lead to impractical shim coils with a very high dissipated power. P T max constraints that are too strict may result in poor shimming performance. In order to balance the trade-off between the dissipated power and the shimming capability, P T max was individually selected for every field map to obtain a dedicated shim coil with the same shimming performance, i.e. the same standard deviation of the residual field, as the unconstrained SH shimming with a given order. In this manuscript, this given order will be referred to as target order and is shown in Table 1.

| Singular value decomposition
For each field map in the training set, a stream function was obtained as described above. All stream functions were used to form a n × m matrix where n indicates the total number of nodes on the mesh and m corresponds to the total number of field maps entering the basis generation algorithm. A singular value decomposition (SVD) was performed on the matrix, and the resulting left-singular vectors (referred to as components) corresponding to a certain number of largest singular values were used to obtain the wire patterns of coil elements of a shim coil array. For each component which can be used as a stream function, field maps per unit current were then calculated by formula (2) and the Biot-Savart law to form a sensitivity matrix. To test the performance of the designed shim array, for each brain field map, the currents flowing through all coil elements were obtained by minimizing the standard deviation of the residual magnetic field over the whole brain.
In the above approach, there is still an open question related to the number of components needed for the coil array. To address this challenge, we gradually increased the number of components used to shim each measured field map until the same or better performance was achieved as the unconstrained SH shimming with a certain order. This type of order is referred to as an achieved order (Table 1). P A max (Table 1) is used to indicate the highest dissipated power of all the used SVD components to obtain a given achieved order SH shimming capability for a field map. The value of P A max may be different for each field map.

| Leave-one-out cross validation
The optimal shim design will be achieved if all measured field maps (12 subjects × 7 positions) in the training set are included in the SVD method. However, it is important to assess the stability of this optimal solution and whether it will perform well in other subjects whose field maps are not included in the design. This property is called the generalization ability of a model in machine learning. 28 In order to test the generalization ability of the resulting shim coil array, a leave-one-out cross validation 20 was used. In the leave-one-out cross validation technique, out of 12 subjects, 11 subjects (11 × 7 maps) were assigned to the training group and 1 was assigned to the validation group. All field maps of the training group were used as target fields in problem (1) and the resulting m = 11 × 7 = 77 stream functions entered the SVD design algorithm from Section 2.3. Thereafter for each field map of the validation group, a minimum number of components was determined for the given achieved order, as described in Section 2.3.
The entire leave-one-out cross validation procedure was repeated 12 times, where each subject was sequentially assigned to the validation group. For example, initially, the first subject constituted the validation group and for each of the seven field maps corresponding to the seven different head positions the number of SVD components from the current training group was defined. Thereafter, the second subject was assigned to the validation group and another seven minimum required numbers of components from the new training group were obtained. Note that in these two steps, the SVD components and correspondingly the shim coil designs are generally different. Finally, 84 minimum required numbers of components were obtained and their statistical parameters (maximum, mean, median, etc.) were considered to assess T A B L E 1 Definition of some useful terms

T n
The total number of required singular vectors (=shim channels) The maximal allowed dissipated power as a constraint parameter in the optimization process of the stream function. This value is chosen for each field map individually.

Target Order
The target order is a chosen SH order which is used to select the maximal dissipated power ( P T max ) as a constraint in the optimization process of the stream function. Strictly speaking P T max is chosen to obtain the same standard deviation in the residual field as that of an unconstrained SH shim system with the same order. (=the target order)

Achieved Order
The achieved order corresponds to the equivalent shimming performance of an SH shimming system with unconstrained power dissipation. For a SVD shim system with a given achieved order, the minimal number of singular vectors (=channels) were chosen to obtain the same SH-shimming performance.
The highest dissipated power of all the SVD channels to obtain a certain achieved order SH shimming capability for a field map.
performance of the proposed SVD algorithm for the given achieved order.

| Shimming capability and coil element layouts of the designed coil array
As described in the above subsection, 84 separate required numbers of components were calculated after the leaveone-out cross validation. The maximum of these 84 minimum required numbers was used as the final total number (T n ( Table 1)) of elements in the coil array design. The SVD algorithm was executed for all 84 measured field maps in the training set, as described in the Section 2.3, and T n components corresponding to the largest singular values were selected. Winding patterns of the channels of the final shim coil array were produced as contours of the respective SVD components. 19 All coil elements were then applied to perform shimming for each field map in the training and test set. The standard deviation of the residual field inhomogeneity after shimming for each map was calculated in order to assess the shimming performance of the resulting coil array.

| Bootstrap approach for testing the stability of the designed coil array
In the above subsections, a shim coil array design method has been proposed using all measured field maps of 12 subjects in the training set which are only a small sample of the population of all human subjects. In the following, we apply an inferential statistical analysis to deduce properties of the population from this sample. In this subsection, we first define a continuous parameter related to the shimming capacity of the designed coil array. Thereafter, the confidence interval of this parameter is calculated and the stability of the designed coil array is analyzed based on that confidence interval. The new parameter is needed because previously used parameters are either discrete integers (e.g. minimum number of components, achieved order, etc.) or contaminated by a bias (e.g. SD of the residual is contaminated by the unshimmable field terms due to local sources 5 ) and therefore are unsuitable for statistical analysis. The parameter is defined as follows: First, for each field map of each subject, shimming using the designed coil array and the target order SH shimming is obtained. The standard deviation of the difference between these two shim fields is then calculated for each field map. Finally, the mean of the standard deviations for all the field maps in the population is specified as the parameter, which is denoted μ σ . Its corresponding statistical estimate for the sample of the 12 subjects is expressed as μ σ . The parameter is used to assess the stability of the designed coil array.
While the calculated parameter μ σ is valid for the given sample of 12 subjects, its distribution over general samples of field inhomogeneities in the brains of all human beings is unknown. In order to obtain an approximate confidence interval of the parameter μ σ using the field maps of the 12 measured subjects, we used a bootstrap approach 21 as described below. First, k subjects, namely a bootstrap sample, are obtained by randomly resampling with replacement from the 12 original subjects. Here, k may not necessarily be equal to 12 and two or more subjects may be repeated. Second, all measured field maps of the k subjects are used in the SVD method from Section 2.3 to generate new T n components. Finally, these new components are applied to shim all field maps in the training set and μ σ using the new components from the bootstrap sample is then calculated. The whole procedure is repeated b times and μ σ for the bth time is denoted by μ b σ . Here, b is set to 1000, which is a typical number in confidence interval approximations with bootstrap. Thus, sorting all the μ b σ in ascending order, the approximated bootstrap percentile 95% confidence interval 21 of the parameter is given by (μ 25 σ , μ 975 σ ) and the significance level α is set to 0.05. The mean of all the μ b σ is indicated by μ σ . The μ σ obtained by the designed coil array and all 84 measured field maps in the original sample is expressed as μ σ , which is used to perform a comparison with the upper bound of the confidence interval in order to measure the stability of the designed coil array.

| RESULTS
This section consists of five subsections. Example coil layouts for individual brain field maps are shown in the first subsection. The results in the next two subsections are obtained using the SVD method described in Section 2.3 and the leaveone-out cross validation technique described in Section 2.4. Coil element layouts of the designed array and the comparison of shimming capability between the conventional SH method and the SVD method with the selected number of shimming channels are presented in the fourth subsection. The results of the stability of the designed array are described in the fifth subsection. Figure 1A-D shows shim coil layouts for various selected subjects in the reference head positions. Here, P T max in problem (1) was selected for each field map to obtain the same shimming performance as the fifth-order SH shimming (target order 5). It can be seen that all coil layouts show some similarities although substantial differences between some layouts, for example as depicted in Figure 1A,C, are visible. Figure 1E-F presents coil layouts for subject 1 after a rotation around the x and y axes, respectively. Here, the designed coil after the head rotation around the z-axis is  Figure 1A), with the corresponding rotation. Contrary to that, the coil layout after the head rotation around the x-axis shows the largest variation from the reference. The variation may be due to the fact that the head rotation around the x-axis leads to a stronger change of susceptibility-induced field distribution in the brain, as discussed in Maclaren et al. 29 Figures 2A-D and 3 show box plots of minimum numbers and maximum dissipated power P A max of the components, which are used to obtain the third to sixth achieved order SH shimming capability for the original field maps in the training set, respectively. Here, all 84 column vectors of the matrix entering the SVD algorithm were calculated separately F I G U R E 2 Box plot of the number of required components (referred to as channel counts) needed to achieve third-order (A), fourth-order (B), fifth-order (C), and sixth-order (D) SH shimming capability for all field maps in the training set. To provide a reference for the plotted channel counts, the numbers of SH functions corresponding to the given order are displayed as black dashed lines, defined as the total numbers of SH functions corresponding to the given order minus four. The box plot of maximum dissipated power P T max [W] (E) for those target orders reveals that P T max grows exponentially with the increase in target order. All those figures also suggest that higher dissipated power P T max can lead to less channel count from the perspective of the median value of channel counts | 1449 JIA et Al.

| SVD method for shim coil design
for different target orders ranging from the third to the eighth order, as shown on the x-axis of Figures 2A-D and 3. Figure 2E shows the box plot of dissipated power P T max selected to solve Problem (1) for the corresponding target orders. In order to assess the number of required components, the total numbers of SH functions corresponding to the different orders (termed SH numbers) are also displayed in black dashed lines in Figure 2. Note that the SH numbers are defined as the total numbers of SH functions corresponding to the given order minus four, since constant and linear SH functions have already been used to pre-process brain field maps. The variations of the maximum and median numbers of components with respect to target orders are individually depicted in red solid lines and red dashed lines with the star symbols in Figure 2. The outliers of the minimum numbers are plotted using plus symbols.
As can be seen in Figure 2, with increasing target order, the median values of the number of required components tend to decrease. Moreover, if a certain target order is equal to or greater than an achieved order, the median values of the number of components are significantly lower than SH functions corresponding to the achieved order. For example, the median values of all minimum numbers of components are only 50%, 52.4%, 40.625%, and 42.2% of the SH numbers if both the target orders and the achieved orders are equal to 3, 4, 5 and 6, respectively. if the target orders are selected to be one higher than the achieved orders, namely 4, 5, 6, and 7, then the median values of the number of components reduce to 16.67%, 19.05%, 15.625%, and 24.44% compared with the corresponding SH numbers. Although the selection of higher target orders may lead to an increase of dissipated power of the used components, as shown in Figure 3, it seems acceptable from an engineering perspective based on previously realized air-cooled coils. For example, the sixth target order and the fifth achieved order can be used in the coil array design since the maximum values of current and maximum F I G U R E 3 Box plot of maximum dissipated power P A max [W] of required components needed to achieve third-order (A), fourth-order (B), fifth-order (C), and sixth-order (D) SH shimming capability for all field maps. The variations of the maximum and median values of these maximum dissipated power with respect to target order are depicted in red solid lines and red dashed lines with the star symbols, respectively. The outliers of the maximum dissipated power are plotted using the plus symbols power dissipation are 9.98 A and 15.77 W, respectively. In this example, each component has been discretized by its contours with 26 steps.
If a target order is lower than the achieved order, the median values of the number of components may be either slightly lower or higher than the SH numbers, as shown in Figure 2B-D, respectively. Furthermore, as shown in Figure 3, the median values of the maximum power dissipation remain within the same range as those in the case where the target order is equal to the achieved order. These results suggest that it is advantageous to use an equal or higher target order than the achieved order when defining the total number of required components.
Although the median values of the numbers of components have a decreasing tendency, the maximum number of components in some cases tend to first decline and then rise with the increase of the target order, as shown in Figure 2A-C. This U-shape variation suggests that using much higher target order than its achieved order could lead to a suboptimal number of required components. In such cases, we can select the transition points, which indicate the points close to the local minimum on the U-shape tendency, to be the number of coil elements of the designed array since the maximum of the power dissipation has an increasing tendency with the target order. For example, as shown in Figure 2B, the target order at the transition point is equal to five and the maximum required number of components is equal to 11, which is the lowest requirement of the number of the coil elements to obtain a shimming capability better than fourth-order SH for all field maps. Figure 2B-D also shows that the maximum number of components can be less than the total number of SH functions if a target order is equal to the achieved order. For example, when both the target order and the achieved order are equal to 5 or 6, the numbers of components for all the field maps in the training set are required to be 24 or 30, respectively. F I G U R E 4 Box plot of channel counts needed to achieve third-order (A), fourth-order (B), fifth-order (C), and sixth-order (D) SH shimming capability in validation groups. To provide a reference for the plotted channel counts, the numbers of SH functions corresponding to the given order are displayed as black dashed lines, defined as the total numbers of SH functions corresponding to the given order minus four However, to obtain the equivalent shimming capability to the SH fifth order (36 components for the full fifth-order SH basis) or sixth order (49 components for the full sixth-order SH basis), for more than 90% of all the field maps, only 19 or 24 components are needed, respectively. Figures 4 and 5 show box plots of minimum numbers and maximum dissipated power of the used components in all validation groups as described in Section 2.4. As can be seen, the median values of the maximum numbers of components also decrease with the increase of a target order. Moreover, if a target order is equal to or greater than an achieved order, the median values of the numbers of components also tend to be lower than the total numbers of the SH functions. For example, the median values of all numbers of components are less than 75%, 71.43%, 75%, and 82% of the corresponding SH numbers when both the target orders and the achieved orders are equal to 3, 4, 5, and 6, respectively. Furthermore, if the target orders are selected to be one order higher than the achieved orders, that is 4, 5, 6, and 7, the median values of the numbers of components reduce to 41.67%, 23.8%, 39.06%, and 55.56% compared to the corresponding SH numbers. These results suggest that using a higher target order at the design stage can enhance the generalization ability of the resulting coil array. Figure 4A-D also shows that the maximum numbers of components in some cases tend to follow a U-shape with the increase of the target order. In these cases, the values at the transition points can be selected as the numbers of coil elements of the designed array. For example, 12 and 24 can be selected as total numbers of coil elements to obtain fourth and fifth achieved order SH shimming capability for all field maps in the training set when the target orders are equal to 5 and 6, respectively. Figure 6A,B shows scatter plots of standard deviations for all field maps in the training set using 12 and 24 components, F I G U R E 5 Box plot of maximum dissipated power P A max [W] of required components needed to achieve third-order (A), fourth-order (B), fifth-order (C), and sixth-order (D) SH shimming capability in validation groups respectively, and compared to high-order SH shimming. As can be seen, for all field maps, shimming using 12 and 24 components outperforms fourth-and fifth-order SH shimming, respectively. Furthermore, for more than 39.29 and 90.48 percent of all field maps, these components achieve an equal or better shimming capability compared to fifth-and sixth-order SH shimming, respectively. Figure 7A-D depicts the first four components of the designed shim array of 12 components. As can be seen, these four components show distinctly different coil layouts in comparison to conventional SH coils. 30 Figure 7E shows a box plot of currents flowing through the 12 components to obtain a fifth achieved order SH shimming capability for all 84 field maps in the training set. As depicted, the median of currents flowing through the first SVD component has a more than three times higher amplitude than those of other components. This result may be related to the fact that the singular value of the first component is greater by a factor of over 2.7 than the singular values of the other components.

| Shimming capability and coil element layouts of the designed coil-array
The performance of the designed coils was then evaluated in a test set consisting of data from 6 further subjects who had not been used during the design phase. Figure 8A presents a box plot of standard deviation (SD) of residual magnetic field inhomogeneity after shimming the field maps in the test set with SH and SVD coil. As can be seen, shimming using the 12 and 24 SVD components performed slightly better than the fourth-and fifth-order SH shimming, respectively. Figure 8B depicts the residual fields on three typical slices of a representative field map in the test set after shimming using the first and fourth SH shapes and SVD coil with 12 elements. As shown, the SVD coil shimming has a similar capability as the fourth-order SH shimming. Figure 9 shows the variation of the 95% confidence interval of μ σ with respect to the number of subjects k in the bootstrap. Here, μ σ is calculated using the designed coil array with 12 coil elements. In order to assess the stability of the designed coil array, μ σ for all field maps in the original sample of 12 subjects is depicted by a dashed black straight line. The upper part of the confidence interval tends to decrease with increasing k. When the number of subjects k is equal to 12, the upper boundary of the confidence interval of μ σ is 20.75 % larger than μ σ , indicating that a different sample of 12 subjects may lead to quite different designs and shimming capabilities. In order to bring the increase of the upper boundary under 5 %, the number of subjects k would need to be around 140.

| DISCUSSION
In this study, we have developed a novel approach for a shim coil array design based on the SVD method with the explicit consideration of the human head anatomy. Instead of directly applying SVD to multiple field maps of human brains, as described in Adalsteinsson et al 16  approach has been recently presented in a design for a gradient coil array with movable imaging volume. 31 The input of the SVD algorithm in our method was composed of individual coil layouts optimized to achieve a given order SH shimming capability for each field map (given a constraint on total power dissipation) with each coil representing a column vector of the matrix. Contours of each resulting component can be used to construct an individual coil element of the shim coil array. Unlike the results of the SVD algorithm field maps in Li et al, 17 the magnetic field generated by each coil element is spatially smooth over the human brain and naturally satisfies Maxwell's equations. Furthermore, the proposed protocol accounts for different possible head positions, which increases the usability of the designed shim array.
In the previous work, 17 leave-one-out cross validation has already been used to evaluate the generalization ability of the SVD modes of the field maps. Here, the cross validation technique has been applied to obtain a suitable total required number of components to balance the trade-off between generalization ability and fabrication challenges of the resulting shim coil array. Furthermore, the power dissipation of each element of a shim coil array can be calculated to provide additional clues for determining the total number of coil elements.
In a less recent publication, 16 it was shown that the first SVD mode combined with up to second-order SH shapes outperforms up to third-order SH functions for the reduction of field inhomogeneities in validation data. This means that six bases (1 SVD mode plus 5 second-order SH shapes) are required if the field maps were measured after only a first-order SH shimming. This number of the bases is comparable to 5, which was the number of coil elements in F I G U R E 8 Box plot (A) of standard deviation of the residual field inhomogeneity after different shimming methods for the data in the test set and the residual fields (B) on three typical slices of a representative field map in the test set after shimming using the first (SH1), fourth SH shapes (SH4), and SVD coil with 12 elements (SVD12TO5) the present work if the target order and the achieved order are equal to 5 and 3, respectively, as shown in Figure 4A. These comparable results may be due to the high degree of similarity of field maps of normal human brains within the general population.
The validation data shown in Figure 4 reveal that using a target order equal to or lower than an achieved order did not allow us to obtain a smaller number of components than that of the corresponding SH functions when the achieved order was higher than 4. This result may be caused by relatively simple resulting stream functions from the individual shim coil design step, according to the method in Section 2.2. The authors 32 have shown that the stream function for a field map has more wire loops for higher SH shimming fidelity. If the target order is too low, the resulting stream function does not capture the complexity of the target field map. After the SVD is performed on the matrix consisting of such stream functions for all field maps, more components are required to obtain a higher achieved order shimming fidelity.
As shown in Figures 2A-C and 4B,C, the maximum number of components tends to rise with the increase of the target order though the median values of the numbers of components have a decreasing tendency. This result may be due to relatively complicated resulting stream functions. If a target order is too high, the resulting stream function for each field map has more local wire loop structures according to Jia et al. 32 Then, the SVD algorithm has difficulties finding a common basis in a set of rapidly oscillating functions. When the resulting SVD components are used for shimming for a specific field map, more components may be required.
As shown in Figures 2 and 4, the total number of used components are typically determined by the maxima or the maximum outliers in the box-plots. One possible reason for the appearance of the outliers is the difference in the magnetic susceptibility distributions for individual human brains. Another possible reason is the different rotations for the same human brain. For example, as shown in Figure 6B, using 24 components did not allow us to achieve a sixth-order SH shimming performance for around 10 percent of the field maps. Even for the same subject, the sixth-order SH shimming capability for field maps corresponds to some but not all head positions. In order to achieve the goal for all field maps, one possible solution is to increase the number of used components to 30, as shown in Figure 2. The relationship between the head rotations and the number of the required components is a topic of an ongoing study at our institution.
It is worth noting that our results are obtained using field maps for whole brain shimming. Although the resulting shim coil array can also be used for reducing the field inhomogeneity within local regions in the brain, it may not be optimal in this case. Moreover, the proposed coil array may not be suitable to perform shimming over other human organs which have different susceptibility distributions. Fortunately, the proposed design method can be applied to create specific shim coil arrays for the local regions or other organs if field maps or the corresponding organs or areas are used in the design.
Stockmann et al 2 has simulated the shimming performance of their integrated RF/B 0 helmet arrays for human brain field maps acquired after second-order global SH shimming. Based on their results, 32-channel shim arrays plus 8-channel face loops have a shimming performance comparable to between the fourth-and fifth-order SH shimming. However, Figure 4C shows that the 24-channel coil elements achieve similar shimming fidelity for all validation data as the fifth-order SH shapes. This means that the total number of coil elements can be reduced by a factor of 1.25 compared to the RF helmet design. This reduction may be due to the fact that two shim coil arrays are mounted on different current-carrying surfaces and our shim coil array was designed specifically to reduce the field inhomogeneity of the human brain.
As can be seen in Figure 6A, some field maps in this study show higher standard deviations than some examples from the literature (e.g. Stockmann et al 12 ). One of the possible reasons of the high standard deviation observed for some field maps in our data set is the large amount of anatomical diversity amongst the subjects enrolled in this study. Another reason can be the selection of the brain region to be shimmed, which in our study comprises the whole brain, including cerebellum and the most part of the brain stem. Furthermore, in our study we have paid a F I G U R E 9 95% confidence interval of μ σ using the bootstrap technique. Here, μ σ is the mean of the standard deviation of the difference between the two resulting shimming fields generated by the designed coil array with 12 coil elements and the target order SH shapes for shimming of the residual field inhomogeneities great deal of attention to accurate segmentation of known regions of strong field inhomogeneities, such as temporal lobes and prefrontal areas.
This study only focuses on global static shimming of the entire human brain. However, dynamic slice-wise shimming of the human brain has attracted increasing interest in the past few years 2,8,10,12,33 since dynamic shimming typically outperforms global shimming. Unlike multi-coils with local loops, 8,10,12 each coil element of the resulting shim coil array has a global distribution of wire loops on the currentcarrying surface. Like SH shim coils, 33 this global distribution might lead to eddy current-induced field perturbations on the B 0 field during dynamic shimming. One possible solution to address this problem is to use a similar technique of pre-emphasis and B 0 compensation as in Juchem et al. 33 There are several challenges associated with the practical implementation of the resulting coil designs: (i) current SVD components are all defined on the same current-carrying surface and (ii) do not account for the practical manufacturing constraints such as a finite wire spacing. Finally, as shown in Figure 7, the topology of the wire patterns is very irregular. The latter issue may be addressed by automated wire layout calculation routines, one such algorithm is being currently developed in our lab. Combined with machined manufacturing technologies, such as etching or water-jetting, the complexity of the underlying winding pattern can be effectively handled. The other two issues may be simultaneously handled by the following procedure: the field maps of all SVD elements can be calculated for a reasonable common region, e.g. an elliptical region encompassing typical brain dimensions, and used as target fields for finding the final stream function values on a set of concentric cylindrical surfaces. During this final design step, the p-norm (‖ ⃗ J‖ p : = ( ∫ Γ � ⃗ J� p dΓ) 1∕p (p > 2)) of a surface current density 34 can be used to control the minimum width of wire tracks, leading to realizable coil element layouts. Following the implementation of these steps the final comparison to SH shimming with a larger number of SH functions will become feasible, as only then it will become possible both to compare the resulting electrical parameters and to assess the practical fabrication challenges.
The manuscript presents a novel shim coil design concept. The current results can be used as a preliminary data to guide further research work. For example, in order to increase the stability of the designed coil array, more field maps from a high number of subjects are needed during the design of a coil array, as shown in Figure 9. In order to reduce the cost of obtaining many field maps, the existing field maps from other projects, such as the human connectome project, 35 will be used in the future. Moreover, the use of a high number of field maps may become a challenge due to increased computational costs. Parallel computations using graphic processing units (GPUs) may be a solution that will be explored in future work.

| CONCLUSION
A novel approach based on the stream function method and the singular value decomposition has been presented and used to design a shim coil array matched to the human brain anatomy. The cross validation technique has been applied for finding an optimal number of elements of the coil array. The results demonstrate the validity of the proposed method to design a potentially physically realizable shim coil array matched to the human brain anatomy, which naturally satisfies the laws of electrodynamics. A bootstrapping technique has been leveraged to estimate the required number of brain field maps for the final coil design.