Robust deep learning estimation of cortical bone porosity from MR T1-weighted images for individualized transcranial focused ultrasound planning

Abstract

Objective: Transcranial focused ultrasound (tFUS) is an emerging neuromodulation approach that has been demonstrated in animals but is difficult to translate to humans because of acoustic attenuation and scattering in the skull. Optimal dose delivery requires subject-specific skull porosity estimates which has traditionally been done using CT. We propose a deep learning (DL) estimation of skull porosity from T1-weighted MRI images which removes the need for radiation-inducing CT scans.

Approach: We evaluate the impact of different DL approaches, including network architecture, input size and dimensionality, multichannel inputs, data augmentation, and loss functions. We also propose back-propagation in the mask (BIM), a method whereby only voxels inside the skull mask contribute to training. We evaluate the robustness of the best model to input image noise and MRI acquisition parameters and propagate porosity estimation errors in thousands of beam propagation scenarios.

Main results: Our best performing model is a cGAN with a ResNet-9 generator with 3D 64×64×64 inputs trained with L1 and L2 losses. The model achieved a mean absolute error of 6.9% in the test set, compared to 9.5% with the pseudo-CT of Izquierdo et al. (38% improvement) and 9.4% with the generic pixel-to-pixel image translation cGAN pix2pix (36% improvement). Acoustic dose distributions in the thalamus were more accurate with our approach than with the pseudo-CT approach of both Burgos et al. and Izquierdo et al, resulting in near-optimal treatment planning and dose estimation at all frequencies compared to CT (reference).

Significance: Our DL approach porosity estimates with ~7% error, is robust to input image noise and MRI acquisition parameters (sequence, coils, field strength) and yields near-optimal treatment planning and dose estimates for both central (thalamus) and lateral brain targets (amygdala) in the 200-1000 kHz frequency range.

Keywords: Deep learning; MRI; image-to-image translation; neural network; porosity; transcranial focused ultrasound.

Fig. 1.

CT, porosity and T1-weighted MRI image sections for 3 of the 15 subjects included in this study. The porosity map is derived from CT as explained in the text. Red arrows point to high-porosity features corresponding to voxels containing primarily bone marrow, which are visible on both the porosity map and the saturated T1-weighted MR image. In other words, the T1-weighted MR image contain the information necessary to extract the porosity information in a patient-specific manner but doing so is not straightforward because of complex MRI contrast around the grey matter, white matter and the dura.

Fig. 2.

Deep learning configurations tested in this work. A: The base DL model used throughout this work is a conditional generative adversarial network (cGAN) with pixel and discriminator losses. All other models are based on this cGAN, with variations emphasized in red. B: cGAN with cropped inputs (size tested: 64×64, 128×128 and 256×256 for 2D, 64×64×64 for 3D). Cropping inputs yields a smaller generative network and increases the size of the training dataset, thus reducing overfitting. C: Use of the pseudo-CT as an additional input channel. D: Auto-context models are obtained by daisy-chaining cGANs together, which may improve the accuracy of the estimation at the last stage. E: Addition of perceptual loss to the backpropagation, in addition to the cGAN pixel and discriminator losses. F: Backpropagation in the mask only includes pixels located inside the skull mask in the backpropagation process. This is done in an effort to focus the degrees-of-freedom of the network on estimation of skull voxels, since porosity outside the mask is known and trivial (=1).

Fig. 3.

Violin plots of the mean absolute error (MAE) in the test set for various DL variations that were found to significantly impact the accuracy of porosity estimates. The bold numbers are the mean MAE for each variation, while the italic numbers above and below indicate the 99^th and 1^st percentiles of the error distribution, respectively. The grey boxes indicate the implementation variant retained in the final model. A: Impact of network architecture and input type. MRI + pCT indicates that both the MRI and the pseudo-CT (Izquierdo et al.) are used as input to the models (two channels). All inputs are 2D with sizes 256×256. B: Impact of input size and dimensionality. The 3D ResNet with 64×64×64 inputs has the best performance, therefore we retain it in the final model. C: Augmenting the data with rotations systematically improved the robustness of the estimation, whereas image flips had minimal impact. Nevertheless, we chose to perform both rotations and flips during training of the final model (comparisons performed for the 2D 256×256 UNet). D: Limiting the backpropagation from pixels in the skull mask improved the error distribution, therefore we retained this strategy in the final model (comparisons performed for the 2D 256×256 Unet).

Fig. 4.

Violin plots of the mean absolute error (MAE) in the test set for various DL variations that did not significantly impact the accuracy of porosity estimates. The bold numbers indicate the mean MAE for each variation, while the italic numbers above and below indicate the 99^th and 1^st percentiles of the error distribution, respectively. The grey boxes indicate the implementation variant retained in the final model. A: Impact of porosity v square-root of the porosity as target estimation metric. Using the square-root of the porosity theoretically boosts small porosity values in the training process but in practice does not affect the training error much. B: All loss functions evaluated in this work yielded similar network estimation error across the test set, therefore we used a simple L1+L2 norm in the final implementation. C: Auto-context modeling did not improve estimation accuracy for our task, therefore we did not retain this technique in the final model (ACM0).

Fig. 5.

Impact of training set size on model performance. A: Validation curves for the baseline model (2D UNet, 256×256 inputs) trained with an increasing number of subjects. There are N=15 subjects in totals, one is set aside for validation while one is set aside for testing. B: Violin plot of the MAE error in the training set, as a function of the number of subjects in the training set. There is no significant performance improvement for N greater than 9, suggesting that using N=13 is well matched to the complexity of the model.

Fig. 6.

Final model performance, compared to pCT-based porosity estimation of Izquierdo et al. and the generic pix2pix GAN.

Fig. 7.

Representative examples of porosity maps estimated with the pCT (Izquierdo et al.), pix2pix, proposed DL approaches and compared to CT (reference). Also shown on the left are whole-FOV and zoom MRI slices, for reference.

Fig. 8.

Performance of the proposed DL approach as a function of noise level in the input MRI image. The moderate, high and very high noise scenarios correspond to grey matter SNR values of 31, 10 and 4, respectively. For each noise level, we show zoomed panels of the noisy MRI inputs and the porosity image outputs. The mean absolute error (MAE) of porosity estimates is equal to 8.4%, 8.4% and 9.3% in the three noise scenarios, indicating that our proposed DL approach is robust to image noise.

Fig. 9.

Porosity estimation using the proposed DL applied to subjects in the open-source Connectome and CERMEP databases, which were not included in the training, validation and testing. T1-weighted images from those databases were acquired on different MRI systems and with different sequence parameters than the inputs used for training, validation and testing. Estimates are reasonable for all subjects, with high porosity values corresponding to marrow-filled pores that are clearly visible on both the MRI and porosity images (green arrows).

Fig. 10.

Scalp maps of the acoustic intensity deposited in the left thalamus of the test subject (arbitrary units), computed using mSOUND at 200 kHz, 500 kHz and 1000 kHz using porosity maps derived from CT (reference), the proposed DL approach and the pseudo-CT methods of Burgos et al. and Izquierdo et al.