Simultaneous imaging of hard and soft biological tissues in a low-field dental MRI scanner

Magnetic Resonance Imaging (MRI) of hard biological tissues is challenging due to the fleeting lifetime and low strength of their response to resonant stimuli, especially at low magnetic fields. Consequently, the impact of MRI on some medical applications, such as dentistry, continues to be limited. Here, we present three-dimensional reconstructions of ex-vivo human teeth, as well as a rabbit head and part of a cow femur, all obtained at a field strength of 260 mT. These images are the first featuring soft and hard tissues simultaneously at sub-Tesla fields, and they have been acquired in a home-made, special-purpose, pre-medical MRI scanner designed with the goal of demonstrating dental imaging at low field settings. We encode spatial information with two pulse sequences: Pointwise-Encoding Time reduction with Radial Acquisition and a new sequence we have called Double Radial Non-Stop Spin Echo, which we find to perform better than the former. For image reconstruction we employ Algebraic Reconstruction Techniques (ART) as well as standard Fourier methods. An analysis of the resulting images shows that ART reconstructions exhibit a higher signal-to-noise ratio with a more homogeneous noise distribution.


Results
"DentMRI -Gen I" (Fig. 1a) is our first-generation MRI dental scanner, designed and built to demonstrate hard tissue imaging techniques in low magnetic field settings. Operating at a field of ≈ 260 mT provided by a permanent magnet, all parts are standard commodities, 3D-printable, or inexpensive to machine. The total cost of all components adds up to ≈ 150 k €, a small fraction of the price of high field scanners (based on super-conducting magnets) previously used for dental imaging 8,12,13,16,21 .
"Gen I" is meant for technology-demonstration purposes rather than for immediate clinical validation. We have therefore restricted our use to ex-vivo image acquisitions of samples accommodated in a cylindrical field of view of height ≈ 100 mm and diameter ≈ 45 mm, even if the gap between magnet poles is significantly larger. We also enclose the samples inside a Faraday cage to isolate our single radio-frequency (rf) coil from the magnetic www.nature.com/scientificreports/ gradient coils and electromagnetic interference noise present in the laboratory, which is otherwise unshielded. Further details can be found in the "Methods" section. In the following, we will first provide images in which hard and soft biological tissues are simultaneously visible, we will then show 3D reconstructions of bare human teeth, and we will end the section with quantitative comparisons between the performance of PETRA and DRaNSSE pulse sequences, and between Fourier and ART-based mathematical reconstructions (See "Methods" section for further details). Table 1 contains the parameters used for all images in this section.
Combined tissue imaging: rabbit head. For initial demonstration purposes, we first show ex-vivo images of a rabbit head, which was soft boiled in tap water to delay tissue deterioration (Fig. 2a). Figure 2c contains selected slices from the full 3D reconstruction employing a PETRA sequence. The field of view is 46 × 54 × 42 mm 3 , with an isotropic resolution of 0.5 mm. We excite the sample with a hard rf pulse of 10 µ s for Table 1. Image acquisition parameters. "NA" stands for "not applicable". www.nature.com/scientificreports/ a flip angle of ≈ 58 degrees, and the repetition time is TR = 10 ms. These parameters were used in two independent acquisitions: one with a short dead time (90 µ s, limited by ring-down in our setup) to read in the combined signal from hard and soft tissues (Fig. 2c top); and one with a long dead time (1 ms) to remove the short-lived contribution of teeth and skull tissues (Fig. 2c middle). Each radial acquisition lasted 710 µ s (4 ms) in the combined (soft) tissue scan, with a total of 4098 (4098) radial spokes, corresponding to an undersampling factor of 7 (7), and 848 (3904) single points to fill the k-space gap, requiring 75 (25) averages and a total scan time of 61 (33) min for the images in Fig. 2c. Both images are reconstructed with ART using = 0.5 and 2 iterations in Eq. (6), and denoised with Block-Matched Filters 22 . Reconstruction times take ≈ 28 min in a low-end GPU (see "Methods" ). Finally, we subtract one of the above scans from the other to produce the bottom image in Fig. 2c, where only hard tissues are highlighted. Even this basic post-processing is enough to identify the upper (A), bottom (B), and inner (C) teeth, as well as the rabbit jaw (D).
Ultra-short T 2 tissue imaging: bare human teeth. The results shown in Fig. 2 demonstrate the capabilities of "DentMRI -Gen I" for imaging of soft and hard tissues simultaneously at sub-Tesla fields. Since human dental structures differ significantly from those in rabbits, below we demonstrate our scanner's performance for a sample consisting of four human teeth embedded into a piece of ham emulating the gum (Fig. 3b). Figure 3a shows 2 dimensional slices obtained from a PETRA scan and applying ART. We acquired the image with a field of view of 40 × 40 × 40 mm 3 and a nominal isotropic resolution of 0.5 mm. To excite the sample we used a hard pulse of 9 µs to produce a flip angle of ≈ 67 degrees. The dead time was set to 100 µ s and the acquisition time was 650 µ s with a repetition time of 15 ms. We acquired a total of 2546 radial spokes, corresponding to an undersampling factor of 8, and 912 single points. To increase the SNR, we averaged over 70 scans with a total scan time of ≈ 65 min. This image is reconstructed with ART using 2 iterations and = 0.5 in Eq. (6), and denoised with Block-Matched Filters 22 . Reconstruction takes ≈ 10 min in a low-end GPU (see "Methods" ). In Fig. 3a we identify the tooth and ham, visible in the photograph in Fig. 3b. The absence of pulp and nerve means there is no bright structure inside the dentin.
All in all, Fig. 3 demonstrates that human teeth can be imaged with high resolution at low magnetic fields, even if our scan times so far are hardly compatible with clinical use. We examine possible approaches to overcome this limitation in the Discussion section.

PETRA and DRaNSSE.
As shown in Figs. 2 and 3, tissue contrast with T * 2 weighting with PETRA is possible by image subtraction. This is a lengthy procedure, however, which can be significantly improved with DRaNSSE (Fig. 7b). Here we compare one approach against the other.
For these tests we use a piece of a cow femur (Fig. 4a), which is mostly composed of only two tissues (bone, with T 2 ≈ 1 ms, and marrow, with T 2 > 50 ms) and therefore facilitates T 2 discrimination and image analysis. In Fig. 4b, images labeled as "PETRA long (short) t d " correspond to an individual PETRA scan with long (short) dead times, where marrow (and bone) appear visible. For DRaNSSE we run a single scan, which can be used to reconstruct marrow (and bone) from the second (first) acquisition in the sequence, corresponding to images with the label "DRaNSSE SE2 (SE1)". Importantly, the total scan times for this study are kept the same, i.e. the sum of both PETRA scan durations is very close to the single DRaNSSE scan time ( ≈ 65 min). Also common to both sequences are: a flip angle of ≈ 90 degrees, a repetition time of 50 ms, a field of view of 46 × 48 × 32 mm 3 , an isotropic voxel resolution of 1 mm, a sampling rate of 24 kHz (8 kHz for long t d PETRA), and a total of 4896 radial spokes in k-space, which www.nature.com/scientificreports/ is fully sampled in this case. For the short (long) t d PETRA acquisition we set the dead time to 85 µ s (1 ms), the radial acquisition time to 915 µ s (2 ms), we fill the center of k-space with 40 (1528) single points, and we average over 7 (7) acquisitions to increase the SNR with a total scan time of 29 (37) min. For DRaNSSE, TE 1 is set to 60 µs and TE 2 to 10 ms, with a common acquisition time of 1 ms. The overall DRaNSSE scan time was 65 min for 16 averages. Images are reconstructed with ART using 7 iterations and = 0.5 in Eq. (6). Reconstructions takes ≈ 35 min in a low-end GPU (see "Methods" ). It is apparent from a qualitative comparison between the raw (unfiltered) image sets IV and VI in Fig. 4b, that DRaNSSE reconstructions feature a higher SNR than with PETRA. Figure 4c shows the SNR along the red dotted lines in Fig. 4b for four different cases, quantitatively reinforcing this observation (further details on these calculations can be found in "Methods" ). Voxels corresponding to bone and marrow tissue both feature a stronger SNR for DRaNSSE when reconstructed with ART (unlike with Fourier methods, which we discuss below). This is consistent with the fact that signal acquisition in DRaNSSE starts when spins are in full phase coherence after the refocusing pulse, while in PETRA the finite dead time means that spins have already started to dephase before data acquisition. Consequently, for SE 1 echo times comparable to the dead time in PETRA, a stronger signal is expected in DRaNSSE, which also makes it possible to use an echo time TE 1 shorter than the shortest dead time in PETRA. Furthermore, the simultaneous acquisition of both images in DRaNSSE allows for more averaging in the same total scan time. In the scans in Fig. 4 we acquired 16 averages with DRaNSSE, versus 7 for both short and long dead times with PETRA. In addition to the number of averages, the sampling rate also influences the noise level. We used 24 kHz for images IV, VI and VII and 8 kHz for image V. All in all, the SNR is a factor ≈ 1.4 higher with DRaNSSE than with PETRA for soft tissues, and ≈ 2.2 for hard tissues. www.nature.com/scientificreports/ A similar comparison between DRaNSSE and PETRA acquisitions of a rabbit head shows that tissue contrast is also enhanced with the former. Figure 5 shows ART reconstruction slices from DRaNSSE (top) and PETRA (bottom) acquisitions. The right (left) column reconstructions show soft (and hard) tissues. Common to both sequences are: a total scan time of ≈ 30 min, a flip angle of ≈ 90 degrees, a repetition time of 50 ms, a field of view of 44 × 52 × 42 mm 3 , an isotropic voxel resolution of 1 mm, a sampling rate of 26 kHz (5.2 kHz for long t d PETRA), and a total of 1426 radial spokes in k-space, corresponding to an undersampling factor of 5. For the short (long) t d PETRA acquisition we set the dead time to 90 µ s (1 ms), the radial acquisition time to 910 µ s (4 ms), we fill the center of k-space with 64 (496) single points, and we average over 12 (9) acquisitions to increase the SNR with a total scan time of 15 (15) min. For DRaNSSE, TE 1 is set to 60 µs and TE 2 to 10 ms, with a common acquisition time of 1 ms. The overall DRaNSSE scan time was 31 min for 26 averages. All images are produced using = 0.3 and 7 ART iterations. Reconstructions take ≈ 12 min in a low-end GPU (see "Methods" ).
The red arrows in Fig. 5 show the upper (A) and inner (C) teeth (see Fig. 2b). The dead time of 1 ms in PETRA (image XI) is not long enough to fully remove the hard tissue signal, making it difficult to identify hard tissues from a comparison with X. One possible solution is to increase the dead time, but this requires longer overall scan times. For instance, if we double the dead time to t d = 2 ms (still significantly shorter than the 10 ms echo time of SE 2 in IX), then we need to increase the number of single points at the center of k-space by a factor 2 3 . This adds up to 3968 single points keeping the rest of the above settings, prolonging the total scan time to ≈ 50 min. A second possibility is to increase the readout window proportionally to the dead time, i.e. going to an acquisition time of 8 ms in our example. However, such a long acquisition time is affected by T * 2 decay, which degrades the sharpness of the image contours. This effect is already visible in XI for an acquisition time of 4 ms. Previously explored alternatives to short-T 2 imaging 23 are the acquisition of a second gradient echo by flipping the gradient polarity 17 , or the reduction of the number of single points by means of a ramped hybrid encoding 24 . However, they both lead to artifacts in our setup due to eddy currents induced during gradient switching. Another strategy is to acquire different sets of radial spokes with different gradient strengths, a technique known as WASPI www.nature.com/scientificreports/ (water-and fat-suppressed solid-state proton projection imaging) 25 . Our attempts to image with WASPI also led to artifacts, in this case due to T * 2 effects (k-space points at the radial boundary between acquisitions are brighter in the strong gradient configuration).
A further advantage of DRaNSSE with respect to PETRA is that the relevant decoherence time constant is T 2 rather than T * 2 , as a result of the induced spin echoes. This means we can wait until hard tissue signals have been strongly suppressed. Consequently, a comparison between images VIII and IX, where SE 2 = 10 ms, shows a starker tissue contrast in regions A and C than between images X and XI.
Finally, a concluding remark regarding DRaNSSE: while it clearly performs better than PETRA in our setup, we occasionally observe reconstruction artifacts originated by stimulated echoes. One possible strategy to overcome this is to tailor the timings so as to not capture them, although it might prove inconvenient. As usual in spin echo sequences, it is therefore relevant to carefully calibrate coherent manipulations in the setup, to minimize the detrimental effect of stimulated echoes.

Fourier transform vs algebraic reconstruction techniques.
Along with the choice of pulse sequence, the mathematical method employed strongly impacts the reconstruction quality. We used ART in most of the reconstructions shown so far, since we find it performs better than traditional Fourier analysis with our radial k-space sampling. Here we look into this topic in more depth.
The results in Fig. 4b allow for a first comparison. One obvious conclusion is that images IV and VI differ significantly in the encircled regions, where Fourier techniques artificially darken the area, whereas the signal intensity with ART is more homogeneous. We only observe this for low SNR conditions in the system, and we do not have an unambiguous explanation to this effect (this is not due to the excitation profile of the hard rf pulse, we do not use regularization terms in the iterative process, and our Encoding Matrix does not account for gradient non-linearities, magnetic field inhomogeneities or dephasing/relaxation effects). A second visible effect is that the contrast between the small lump of soft tissue and the external surface of the bone where it is attached (white arrows), is higher with ART than with Fourier transformation. Figure 6 shows noise and SNR maps for ART and Fourier reconstructions of the cow femur ( Fig. 6a with DRaNSSE and Fig. 6b with PETRA) and rabbit images (Fig. 6c,d both with PETRA and 1 mm and 0.5 mm, respectively), following procedures described in the "Methods" . We observe significant differences between the noise pattern obtained from Fourier reconstructions respect to ART. The former present a highly non-stationary noise for all images (i.e. the variance of noise σ (x) depends on position). This is reflected by the coefficient variation (CV) of noise, as well as the average noise level, both shown in Table 2. We also observe for Fourier reconstructions that the maximum noise level happens at the gradient isocenter position, which is coherent with the radial sampling scheme due to the high sampling density close to the center of k-space. For ART, the matrix size strongly affects the level of non-stationarity, where the noise variation coefficient is 8 times smaller when compared with Fourier reconstruction for Fig. 6d, with double matrix size. All in all, we observe that, while PETRA images reconstructed with ART present a more stationary noise than Fourier reconstructions, the noise obtained with DRaNSSE is highly non-stationary even when reconstructed with ART.
In addition to the noise performance, there are also important differences in the SNR. Figures 6a,b show that ART reconstructions for both DRaNSSE and PETRA produce higher SNR than Fourier reconstructions in the hard tissue. Generally, ART produces SNR values that, in the worst case, are similar to those from Fourier reconstructions (see Fig. 6c,d and Table 2). Moreover, note that ART yields better results than Fourier transforms both when k-space is fully sampled (Fig. 6a,b) and when it is undersampled (Fig. 6c,d). Finally, since the results with ART depend on the reconstruction parameters and n it , and there is therefore room for potential improvement with respect to the results we present here, we conclude that ART is better suited for image reconstruction than Fourier analysis with our settings.

Discussion
In summary: we have demonstrated the capability of our new low-cost "DentMRI -Gen I" scanner to simultaneously image hard and soft biological tissues; we have devised a new pulse sequence (DRaNSSE) that, compared to standard short T * 2 sequences such as PETRA, yields higher SNR images and enhanced tissue contrast; and we have shown that iterative techniques (ART) outperform traditional Fourier methods in all quantified metrics, except for the computational time required for the reconstruction. Table 2. Estimated noise and SNR parameters corresponding to images in Fig. 6. SNR in soft and hard tissue on Fig. 6c,d correspond to tongue and inner tooth, respectively. CV stands for coefficient of variation (the standard deviation divided by the mean) and . is the average operator.  www.nature.com/scientificreports/ Dental imaging at low field strengths can be advantageous for reasons beyond cost reduction 14 . For example, magnetic susceptibility effects, which lead to artifacts in the presence of metallic fillings or implants, scale linearly with field strength. Besides, T 1 times are shorter for weaker fields (which allows for faster repetition rates), T * 2 coherence times are longer (somewhat mitigating the SNR loss due to the lower field strength), and specific absorption rate (SAR) limitations scale with the square of the field strength (enabling long echo trains for soft tissue imaging). Also off-resonant effects due to chemical shifts in the Larmor frequencies of different tissues are suppressed (tens of Hertz at 260 mT). Despite these advantages, we are not aware of previous dental MRI systems operated at sub-Tesla fields.
Nevertheless, a prospective "Gen II" scanner must pose solutions to significant remaining challenges if it is to be compatible with clinical conditions. First and foremost, scan times must be reduced (see last column in Table 1). For example, the acquisition times for the human teeth images shown are intolerably long. We therefore plan three major upgrades for "Gen II": i) we will make use of balanced steady-state free precession protocols, which yield optimal SNR values and can be safely used at low magnetic fields 15,26 ; ii) we will use quantum dynamical decoupling techniques such as WAHUHA 27 or CHASE 28 sequences, which can prolong the lifetime of the magnetic resonance signal of hard tissues; and iii) we will perform dual species MRI on protons and 31 P nuclei, since the latter are more abundant in hard biological tissues and they provide complementary information 29 . Additionally, we are starting to investigate the possibility of slice selection with zero-echo time sequences for fast 2D imaging 30 . Aside from shortening scan times, a reduction in the overall footprint and weight would ease scanner siting. The heaviest (and most expensive) component is the permanent magnet ( ≈940 kg and ≈40 k€), which on top imposes the use of a sizable mechanical support structure. For this reason, we will place the new www.nature.com/scientificreports/ magnet with the yoke at the bottom and the poles pointing upwards (forming a "U") so patients can lie on a flat bed. This reduces greatly the complexity and footprint of the required mechanical structure. Another important aspect is the radio-frequency system, which will differ significantly from the "Gen I" configuration. In "Gen II", we will adapt the size and geometry to the human physiognomy, where neck and shoulders place strong boundary conditions, and we will count on dedicated transmit and receive channels, where the latter can be miniaturized and placed in direct contact with the subject's face or teeth increasing significantly the SNR 21 .
Besides dental imaging, a low cost MRI scanner capable of detecting hard tissues and solids may find application in different scenarios, including but not restricted to: head and extremity imaging 14,31 ; the food industry (e.g. in inspection and selection tasks, see 32 ); or rubber degradation control, which is relevant in, for example, the industries of mining or transport 33 .

Methods
"DentMRI -Gen I" Scanner. Main magnet. In MRI systems the quantization axis, spin resonance (Larmor) frequency and, to a large extent, the signal-to-noise ratio of the detected signals are typically determined by a strong homogeneous magnetic field 1 . Here we employ a "C"-shaped permanent NdFeB magnet (Sabr Enterprises LLC, Fig. 1a), providing a main evolution field of ≈260 mT. The field is shimmed down to an homogeneity of ≈20 parts-per-million (ppm) over a spherical region of 150 mm in diameter. The magnet is the heaviest component of the scanner with a mass of around 940 kg for a distance between magnet poles of ≈220 mm. In order to distribute its weight over a ≈5 m 2 floor surface (compliant with typical architectural regulations), we designed and constructed a support structure consisting mostly of profiles manufactured out of Aluminum 6063. Although this was specified to be non-magnetic, it induces field-strength inhomogeneities at the 4000 ppm level. These are mostly linear and we regularly shim down to ≈4 ppm during operation with the magnetic gradient system.
Magnetic gradient system. Linear magnetic field inhomogeneities (magnetic gradients) spatially encode information in MRI setups 1 . Intense gradient fields lead to efficient encoding, i.e. strong resolving power. "Gen I" is equipped with a gradient system capable of reaching strengths > 400 mT/m along any spatial direction, sufficient for sub-mm resolution. The gradient system consists of three pairs of planar coils. Each pair produces a magnetic field pointing in the same direction as the main field, and a gradient in magnitude along a Cartesian axis in the laboratory frame of reference. Each of the z (x, y) coils is formed by one (two) lobes (see Fig. 1a). We fabricated the coils with hollow-tube Oxygen-Free-High-Conductivity (OFHC) copper, with low ohmic losses ( < 50 m� per pair) and the possibility of heat removal by running cooling water through the inner conduct. We can stably operate at 100% duty cycle with a total power consumption < 8 kW. In order to pulse the magnetic gradient fields in short times of 50-100µ s, we wound each coil into only a few loops for a total pair inductance < 100 µH in all three cases. The size and disposition of the coils ensures deviations < 10 % from perfect linearity over a spherical region of 10 cm in diameter. The space between loops is filled with translux D150 epoxy to avoid short circuits and improve mechanical stability. We drive the coils with bipolar amplifiers from International Electric Co. (GPA-400-750), which can ramp from 0 to ±400 A in 100 µ s with our loads. Radio-frequency system. In MRI, rf electronics are required to excite the sample spins and detect their response for amplification, digitization and, ultimately, data processing as required for image reconstruction 1 . We employ a single coil for both sample excitation and signal detection. The solenoid, of length ≈100 mm, ≈52 mm in diameter and with 20 windings with copper wire of ≈0.4 mm diameter, is shown in Fig. 1a. The solenoid includes a gap capacitor to homogenize the current along the coil wire, resonant at ω c ≈ 2π·11 MHz with a quality factor Q ≈ 32 . We keep the Q intentionally low to achieve excitation and detection bandwidths > 1 MHz, compatible with high spatial resolution images. The rf electronics for fine tuning and impedance matching of the resonant circuit are placed next to the solenoid (Fig. 1b). Both the coil and matching electronics are enclosed in conducting Faraday boxes to avoid the deleterious effect of electromagnetic noise in the laboratory and coupling to the gradient coils. We 3D-printed the housing structure out of polylactic acid for mechanical stability (Fig. 1a).
In transmission (Tx) mode, we drive the solenoid from an rf power amplifier (RFPA-4/11-2000 from Barthel HF-Technik GmbH) fed by a direct digital synthesizer on our field-programmable-gate-array (FPGA) board (RadioProcessor-G from SpinCore Technologies Inc.). The RFPA output is low-pass filtered and sent to the solenoid coil (probe) for sample excitation. This scheme is shown in Fig. 1b. A TxRx switch (Barthel HF-Technik GmbH, dead time < 5 µs at 11 MHz) commutes the system operation from transmission to reception (Rx). A low noise amplifier (Barthel HF-Technik GmbH, gain 39 dB and noise factor 1.0 dB) amplifies the analog signal induced by the precessing protons on the coil. After low-pass filtering, we digitize and digitally down-convert and filter the signal directly on the FPGA board.
Control system and graphical user interface. The FPGA board is the main component in the experimental control system. This is plugged into a Peripheral Component Interconnect (PCI) slot on the motherboard of a control computer, and allows us to: i) generate low power rf pulses which are amplified for coherent spin excitation; ii) generate three independent low-frequency and low-power outputs which are amplified and fed to the gradient coils for spatial information encoding; iii) read in, digitize, down-convert and filter the MR signal emitted by the sample; iv) execute all of the above operations synchronously; and v) communicate bidirectionally with the control computer.
We have programmed our own Graphical User Interface (GUI) in Matlab, where we design pulse sequences, set the individual pulse parameters independently, configure data acquisition and filtering settings, and visualize Pulse sequences. MRI pulse sequences are designed to manipulate the sample magnetization: resonant rf pulses lead to coherent rotations, and magnetic-field-gradient (or simply gradient) pulses encode spatial information. The short lifetime of hard tissue signals imposes the use of special-purpose pulse sequences. These typically force a high-intensity constructive interference signal (echo) over very short timescales, or even simultaneously excite the sample and detect its response 3 . Ultra-short echo time (UTE 34 ) and zero-echo time (ZTE 16 ) pulse sequences are prominent examples of the former approach, and SWeep Imaging with Fourier Transformation (SWIFT 13 ) of the latter.
ZTE and SWIFT are better suited than UTE for ultra-short T 2 tissues such as those present in teeth. If the scanner hardware allows for rf excitations with sufficient bandwidth, ZTE is the most convenient choice 35 . Since the "DentMRI -Gen I" scanner has a rather small field of view and we have a fast 2 kW rf power amplifier, we acquired all images above with two variations of standard ZTE pulse sequences (Fig. 7): Pointwise Encoding Time Reduction with Radial Acquisition (PETRA 17 ); and a sequence we have devised to overcome specific k-space (spatial frequency space) coverage and contrast limitations, Double Radial Non-Stop Spin Echo (DRaNSSE).
PETRA . In PETRA (Fig. 7a), a hard rf pulse homogeneously excites the sample after the onset of the encoding gradient fields. Data acquisition starts next, after a dead time usually set by the TxRx switch and the rf coil ring-down. Every acquisition follows a radial direction (spoke) in a 3-dimensional k-space. The center of k-space is not sampled in this way due to the finite dead time, so it is filled in a pointwise manner on a Cartesian grid. Once all radial spokes have been sampled, object reconstruction can be carried out either by inverse discrete Fourier transformation (which requires regridding and interpolation operations), or with Algebraic Reconstruction Techniques (ART [18][19][20]. PETRA suffers from two main constraints. On the one hand, scan times can become exceedingly long when the pointwise-filled gap is large. This is accentuated at low fields, due to the slow ring-down of the Tx coil, which happens with a time constant τ = 2Q/ω c . On the other hand, it is not possible to have T 2 contrast with PETRA, since the echo time (TE) is defined to be zero. Our solution to both problems is our new sequence DRaNSSE (Fig. 7b). When we do acquire with PETRA, we can subtract two independent scans with different parameters for T 2 contrast, following a previously existing scheme with UTE pulse sequences 36 . In one acquisition dead times are short, so hard tissues still appear bright, while in the other we use intentionally long dead times, to ensure hard-tissue contributions fade away before readout. www.nature.com/scientificreports/ DRaNSSE. In DRaNSSE, after the initial rf pulse, two refocusing π-pulses produce two different spin echoes. The first one (SE 1 ) is induced at an echo time (TE 1 ) as short as possible to include contributions from both hard and soft tissues. The second echo (SE 2 ) occurs at a later time (TE 2 ), when the short-lived signal from hard tissues has already faded away. A first acquisition takes place between SE 1 and the second refocusing pulse, corresponding to radial spokes in k-space from k = 0 to the maximum sampled value ( k max ); the second acquisition starts after the second π-pulse and ends at TE 2 , sampling from −k max to k = 0 . The second acquisition could go beyond SE 2 , allowing for higher quality imaging of soft tissues. However, the occurrence of multiple unwanted stimulated echoes due to small experimental imperfections (STE 1 -STE 3 in Fig. 7b) can degrade image quality. STE 1 appears a time τ 1 = TE 1 /2 after the second π-pulse, while STE 2 and STE 3 occur at times τ 1 and 2τ 1 after SE 2 , respectively. Since τ 1 in our sequences is much shorter than τ 2 (the time between the start of the first acquisition and the center of the second refocusing pulse), we find it convenient to acquire the signal before SE 2 to avoid reconstruction artifacts. Note that DRaNSSE is similar to radial multi-echo spin-echo with an echo train length of two 37 , but with half-echo acquisition to minimize TE 1 and avoid contamination from stimulated echoes.
Fourier and algebraic reconstruction techniques. Spatial encoding in MRI relies on inhomogeneous magnetic fields, which provide a Larmor or spin-precession frequency ( ω L ) dependent on the position of the nuclei in the Region of Interest (RoI). Mathematically, this can be expressed as where γ is the gyromagnetic ratio ( ≈ 2π · 42 MHz/T for protons) and B is the magnetic field at position r at time t. As the pulse sequence evolves, the phase acquired by the spins depends on their position: During their precession, spins induce a time-varying signal with the interference of all spins on a nearby detector: where ρ(r) is the spin density distribution in the RoI. This signal is then digitized during a readout or acquisition window, and we apply one of the following mathematical tools to reconstruct an image.
Discrete Fourier transform reconstruction. Typically, scanners make use of linear gradient fields, in the presence of which the Larmor frequency varies also linearly with position. In this scenario, an inverse Fourier Transformation (FT) of s(t) suffices to reconstruct ρ(r) , since the integral in Eq. (2) after down-mixing is trivial and Eq. (3) becomes Here we have assumed, without loss of generality, that a gradient of strength G z points along the z-axis. Since the detected signal is discretized in time by an analog-to-digital converter, fast Fourier transform protocols can be used to map a discrete k-space onto a discrete reconstruction in real space 1 . Discrete Fourier transforms are computationally efficient and applicable in many useful scenarios, but this simple approach may yield suboptimal results when the sampling data does not adjust to a Cartesian grid in k-space. This occurs, for instance, if gradients are time dependent, discretization times are not homogeneously distributed or, as in the present work, k-space is sampled radially or following curved trajectories. In such cases, the acquired data must be pre-processed (e.g. with regridding, density compensation or interpolation operations) prior to Fourier transformation 38 . To perform Fourier reconstruction we first interpolate the k-space data to a Cartesian grid and then apply a Fast Fourier Transform protocol to the interpolated data.
Reconstruction based on an encoding matrix. An alternative to pre-processing and Fourier operations in non-Cartesian k-space sampling is to build a linear forward model for the system response to the applied pulse sequence, define a cost function for the reconstruction, optionally add regularization terms to penalize unrealistic results, and solve a linear inversion problem 39 .
Once the time-dependent signal resulting from the interference of the precessing nuclei has been recorded and discretized, s(t) becomes a vector S of length equal to the number of time steps n t , ρ(r) becomes a vector ρ of length equal to the number of voxels n v , and exp −i�(r, t) becomes the Encoding Matrix with n t rows and n v columns. Equation (3) thus changes to and ρ can be obtained by direct inversion of the Encoding Matrix as −1 S , or by any other means of solving the system of linear equations, e.g. by iterative algorithms such as Algebraic Reconstruction Techniques (ART) [18][19][20] . ART estimates ρ(r) based on the recursive equation   www.nature.com/scientificreports/ where is a control parameter, S i is the i th component in vector S , i is the i th row in the encoding matrix , and ρ 0 can be set to zero. The estimated solution ρ n is updated n t · n it times through Eq. (6), where n it stands for the overall number of ART iterations.
Although iterative methods such as ART are computationally slow compared to discrete Fourier transforms, we find in this work that they can vastly outperform Fourier reconstruction in other relevant metrics.
Estimation of the SNR. We assume that the final reconstruction signals in image domain follow a non-stationary Rician distribution 40 , i.e., a Rician distribution where the noise parameter σ becomes position dependent: σ (x) . This assumption is based on two facts: 1. Gaussianity and Rician: ART and FT reconstructions both produce signals with additive Gaussian Noise in each of the points. The magnitude of that Gaussian data produces Rician data if the real and imaginary parts of the former are independent and have the same σ parameter. Although this need not be strictly the case with our reconstruction processes, we show below that the Rician approach yields valid results. 2. Non-stationarity: ART and FT both perform a linear reconstruction of complex Gaussian data where different weights are applied to the original values. As an effect of the weights and the correlations introduced by the reconstruction, the final reconstructed data can show different noise properties in each voxel, i.e., every voxel can show a different value of σ (x) . In our images, as is usually the case, σ (x) varies slowly with x and can be considered a low frequency signal 40,41 .
In order to test the first assumption (the Rician nature of the data), we consider Fig. 6d, where the background noise is to a good approximation stationary and a single value of σ characterizes the distribution. The expected signal strength in the background is zero, and the Rician Probability Distribution Function (PDF) simplifies to a Rayleigh distribution. From the fits in Fig. 8 we conclude that the signal in the background is Rayleighdistributed, which is compatible with our Rician assumption.
Thus, in order to estimate the SNR of the final volume we assume a non-stationary Rician model. We define the estimated SNR as where S(x) and σ (x) are, respectively, an estimation of the original signal (in absence of noise) and and estimation of the variance of the noise in each point of the image.
In order to estimate the noise, we use a 3D extension of the homomorphic scheme proposed in Ref. 41 . The rough estimation of the original signal is provided by the second order moment of a Rician distribution. If M(x) is a Rician signal, then 40 We estimate the expectation E{M 2 (x)} by a local mean, calculated as the convolution of the original signal with a 3 × 3 × 3 average kernel, h(x) . Hence: and the estimated SNR becomes: , E{M 2 (x)} = S 2 (x) + 2 · σ 2 (x).
S(x) = max{M 2 (x) * h(x) − 2 · σ 2 (x), 0}, Figure 8. Fit of the background data of image I in Fig. 2 to a Rayleigh distribution. A fit to a Normal distribution is shown for comparison.