This is a plain English summary of the following open-access paper:

Marty, P., Boehm, C., van Driel, M., & Fichtner, A. (2024). Transcranial ultrasound modeling using the spectral-element method. The Journal of the Acoustical Society of America, 156(6), 3674–3693. https://doi.org/10.1121/10.0034474

This article assumes that you already have some background in a few wave modeling related topics. If you’re new to wave physics, maybe check out these other explainers first:

• Wave simulations
• Finite-element meshing

TL;DR

This paper provides a framework for modeling the propagation of ultrasound waves through the skull using a technique called the spectral-element method. We show that this approach is extremely computationally efficient and is capable of dealing with highly complex geometries.

Research Questions

There are two research questions that this paper answers:

1. How can the propagation of ultrasound waves be modeled for applications involving the skull using the spectral-element method?
2. How can attenuation (damping of the ultrasound waves) be incorporated into this type of modeling if the domain contains both fluid (soft tissue) and solid (bone) regions?

Motivation and Purpose

Transcranial ultrasound involves modeling the propagation of ultrasound waves through the intact skull. Trying to propagate the ultrasound waves through the skull introduces a number of significant challenges; I’ve detailed these in this previous blog post. Our goal is to develop wave modeling methods that work well for dealing with these complexities.

While imaging the brain is one of the major areas that we’re interested in applying these techniques for, this paper primarily demonstrates these approaches for focused ultrasound applications. Focused ultrasound involves using special types of ultrasound sensors that allow for energy to be concentrated at a specific point in the tissue and can be used to non-invasively treat a variety of different brain-related conditions.

The key benefit of focused ultrasound compared to other existing surgical interventions is that the ultrasound waves can be delivered to the interior brain tissue through the intact skull. This means that it is possible to treat a variety of brain-related ailments that would otherwise require one to undergo invasive open-skull surgical procedures.

A Bit of Background

Here we’ll provide a very brief overview of:

• How waves propagate in transcranial ultrasound.
• How spectral-elements work.
• What attenuation is.

What Do We Actually Model?

Ultrasound waves are just very high-frequency sound waves that propagate through the tissue of the patient. I give a number of examples in this blog post, but here’s one example of what that looks like for the human brain:

Here we can nicely visualize how the ultrasound waves are strongly scattered and distorted by the presence of the skull.

As I mentioned before, we are using special ultrasound sensors here that are designed to deliver considerably more energy into the patient’s tissue when compared to traditional ultrasound devices. In fact, most “normal” ultrasound sensors (think of something like the handheld devices you may have seen in the hospital before) would break if you tried using them for this type of application because they’re just not designed to output this much energy.

What we want to do in focused ultrasound simulations is calculate what the maximum pressure will be that this wavefield produces within the tissue of the patient. This is important for two reasons, namely that the amplitudes at the focus need to be:

1. Sufficiently high to treat the brain-related ailment at hand.
2. Sufficiently low to ensure that the pressure is below some pre-determined safety threshold.

The Spectral-Element Method

The spectral-element method is a technique for modeling the propagation of sound waves. This approach was initially developed in geophysics for modeling domains that are very geometrically complex. Since the spectral-element method is a special type of finite-element method, it allows for the use of meshes that are constructed to resolve some specific boundaries.

This ability to use these types of meshes is one of the key features that makes the spectral-element method so attractive for transcranial ultrasound. Most existing tools that are currently used in transcranial ultrasound are discretized onto rectilinear grids, which means that they split up the domain into a regular grid and then model the propagation of waves on that. While this is a very convenient way of discretizing the model, it has the tendency to introduce so-called staircasing artifacts, unless a very fine discretization is used:

If you want to learn more about staircasing, check out this previous blog post.

Like with other finite-element methods, we can solve the wave equation using a series of matrix operations that are applied over this interconnected network of elements. Generally, the most costly steps for solving this type of system in conventional finite-elements is performing a matrix inversion for the so-called mass matrix.

One of the key benefits of the spectral-element method in wave propagation problems is that this mass matrix is diagonal by construction, meaning that this matrix inversion can be performed (basically) for free. This leads to huge computational savings since one does not need to solve a linear system for each matrix inversion.

In order to preserve the computationally attractive properties of the spectral-element method, we want to maintain this feature. That is, if we make modifications to how we solve the wave equation, we need to make sure that the mass matrix remains diagonal.

Attenuation and Fluid-Solid Coupling

Attenuation is the process by which sound waves get absorbed as they propagate through a certain material. You can think of this like how foam sound panels in a music studio work — the sound panels absorb and dampen the sound waves.

It turns out that the skull is highly attenuative, meaning that it absorbs a lot of the sound that we try propagating through it; this has been well known for several decades. This is a great evolutionary adaptation because it helps to protect our brains from the outside world, but it also makes it much more challenging to deliver ultrasound waves to the interior brain tissue.

Incorporating this type of attenuation into a wave modeling solver in an efficient way is challenging. I won’t go into the details of this here (the underlying math is a bit dense), but the paper goes through how this can be implemented efficiently within the spectral-element method.

The important thing: we’re able to incorporate these attenuating effects in both fluids (soft tissue) and solids (skull) while preserving the computational properties that make the spectral-element method so powerful.

Staircasing Artifact Mitigation

The mitigation of staircasing artifacts by constructing appropriate finite-element meshes is imperative. In order to highlight the importance of mitigating these artifacts in the spectral-element method using boundary conforming meshes, a sensitivity analysis is performed where the maximum pressure distributions are compared for a varying mesh resolutions (number of elements per wavelength).

Two mesh types are compared here, namely:

1. Rectilinear meshes
2. Boundary conforming cubed sphere meshes

Rectilinear meshes are fairly trivial to construct given that they consist of a network of interconnected cuboid-like elements. It is additionally possible to refine the elements along the skull boundary in order to improve how closely we are able to resolve these interfaces. In order to preserve the computationally attractive properties of the spectral-element method, we use local refinements that allow for the conforming nature of the hexahedral mesh to be preserved. We can see an example of what this looks like below:

We want to avoid having small elements in our mesh because we’ll end up taking a hit in the global time step of our wave simulation, which results in the compute cost going up; I discuss this in detail in this blog post. From a computational standpoint, having a boundary conforming mesh is ideal, but constructing this type of mesh is considerably more involved (especially in 3D) when compared to its rectilinear counterparts.

Proof-of-Concept: the Shepp-Logan Phantom

We ran a series of simulations using both these rectilinear and boundary conforming (cubed sphere) meshes to see how greatly these pressure distributions changed between simulations. For this series of tests, we used an idealized skull model known as the Shepp-Logan phantom (from Shepp & Logan, 1974), which consists of an ellipsoidal skull that contains a number of ellipsoidal soft tissue inclusions. While a relatively simple model, it does a pretty good job at highlighting differences introduced by different discretizations.

Here’s a comparison of a few of the maximum pressure distributions within the Shepp-Logan phantom using a focused ultrasound transducer:

Since the cubed sphere meshes, such as in panel (d), are boundary conforming, we are using the cubed sphere with a mesh resolution of 4.0 elements per wavelength (EPW) as the reference solution. The simulation result at 4.0 EPW for the cubed sphere is visually indistinguishable from the one at 2.0 EPW shown in the figure above.

We can quite clearly see the presence of staircasing artifacts along the skull interface within the cases of the rectilinear meshes at 2.0 and 4.0 EPW. While the refined rectilinear grid in panel (c) does a better job at mitigating these staircasing artifacts, it introduces some significant penalties in terms of compute cost, which we’ll discuss below.

We perform a more comprehensive analysis of the different mesh configurations in the paper, but the important take-away is that the boundary conforming cubed sphere meshes achieve excellent results while remaining computationally efficient.

Performance on HPC Systems

All of the above simulations were run on the Alps supercomputing system from the Swiss National Supercomputing Centre (CSCS). At the time of writing, the Alps supercomputer is among the most powerful supercomputing systems in the world and it serves as an excellent platform for testing this spectral-element approach.

The computations that are performed for each wave simulation are distributed across multiple graphics processing units (GPUs) so that they can be computed very efficiently. We decided to allocate around 500000 elements per GPU since we found that this led to a good balance between compute speed and while minimizing inter-GPU communication. Since each compute node on Alps has 4 GPUs, we always rounded the number of GPUs to the nearest integer multiple of 4.

Rectilinear Meshes

Cubed Sphere Meshes

In general, we’ve found that somewhere between 2.0 and 3.0 elements per wavelength gives a good compromise between compute cost and simulation accuracy with fourth-order boundary-conforming hexahedral meshes. It’s also excellent to see that the compute times associated with performing these wave simulations scale very well and that, given a reasonably sized mesh, they can be computed in a matter of minutes.

You may be wondering why the simulation with the refined mesh took so long (over half an hour!) when compared to the other simulations. The reason for this is that the small elements within the areas where we apply these local refinements requires us to use a very small time step to ensure that our simulation remains stable.

So why not just throw more GPUs at the problem to make it go faster? The limitation that we quickly run into with this type of refined mesh is that we’re able to parallelize in space (i.e. take the mesh, split it up into smaller chunks, and send each chunk to one GPU), but the same is not true in time. We are forced to solve this problem sequentially (i.e. how the wavefield looks at the current time step depends on the wavefield from the previous time step), so using more compute resources does not lead to an appreciable speed-up in this case.

The moral of the story: we may be able to get satisfactory results using local mesh refinements, but we’ll take a huge hit in the time step and, thus, a significant penalty in terms of compute cost. This means we should use the boundary conforming meshes in spectral-elements for these types of transcranial applications whenever possible.

Validation Against Existing Models

It’s all well and good to set up some kind of numerical method, but we wanted to validate that what we’re modeling actually gives us physically meaningful results. To check our numerical results, we compared our forward modeling approach to that obtained using the axisymmetric implementation of the k-Wave code.

k-Wave is arguably the most widely used wave solvers in medical ultrasound. It solves the same underlying physical equations that this spectral-element implementation does, but it uses a very different choice of numerical method (known as the k-space pseudospectral method). We chose this tool for validating our results for a few reasons:

1. Since it is so widely used in medical ultrasound, its results have been compared against laboratory measurements within a number of different studies.
2. The axisymmetric implementation we’re considering for this comparison allows for an exceptionally fine discretization to be achieved.

Here we used the same setup as benchmark 6 from the paper by Aubry et al., 2022:

Identical simulation setups are used between k-Wave and the spectral-element code Salvus in order to validate that the solutions computed using the spectral-element method are consistent with other tools commonly used in the medical ultrasound community. Similar to the Shepp-Logan setup, we analyze the maximum pressure distribution throughout the domain, as can be seen below:

We can see that there is excellent agreement between the two solutions. This gives confidence that the solutions computed using the spectral-element method are able to achieve physically meaningful results.

Demonstration on a Complete Skull Model

Both of the previous two numerical examples consisted of relatively simple geometries. However, the degree of flexibility with which the hexahedral finite-element meshes can be constructed is one of the most significant practical advantages of using the spectral-element method for these types of transcranial applications.

For this final example, we consider an adaptation of the multimodal imaging-based detailed anatomical model (MIDA) from Iacono et al., 2015. This numerical phantom consists of the individually segmented tissues in the human head and serves as an excellent model for demonstrating various forward and inverse modeling algorithms.

We first create an unstructured conforming hexahedral mesh using a combination of Blender and Coreform Cubit. We use a similar process to that which we previously published in Marty et al., 2022, where we first construct a high quality surface representation of the soft tissue-skull interfaces using Blender and then create the hexahedral mesh using the meshing from volume fractions algorithm in Coreform Cubit. Here’s the resulting hexahedral mesh that we then use for the subsequent wave simulations:

The material properties that are used for performing the wave simulation are interpolated onto the hexahedral mesh. While this mesh “only” has about 8.5 million elements, the high order nature of spectral-elements means that we store about 1.1 billion nodal values for each of the unique material properties in the mesh (of which we have five in this case).

We wanted to compare what the maximum pressure distributions looked like for cases with and without attenuation. Here’s the result:

We can see that while the overall shape between the two pressure distributions are similar, there are two significant differences:

1. The presence of significantly more scattered arrivals when attenuation is neglected.
2. A considerable decrease (approximately 20%) in the maximum focal pressure when attenuation is included.

It is, however, worth noting that the overall focal shape and volume remains consistent when attenuation is either included or neglected. Here’s a visualization showing the focal volume at -6 dB:

This demonstrates that while the peak pressure and the amount of scattering changes considerably when either including or neglecting the effects of attenuation, the shape of the focus itself remains relatively unchanged. There are a few important implications that we can draw from this:

1. For focused ultrasound, the position and shape of the focus can be approximated quite effectively using a lossless model. However, the peak pressure may be overestimated.
2. For transcranial ultrasound imaging, it is important to include these attenuating effects to be able to correctly fit the lower amplitude arrivals within the ultrasound signals.

Conclusions

We have shown the efficacy and flexibility of the spectral-element method for modeling the propagation of sound waves in transcranial ultrasound. This technique, which is widely used in the field of seismology, is a highly efficient numerical tool for solving the wave equation and excels at dealing with highly complex geometries. The coupled attenuative formulation for spectral-elements that we propose in this manuscript allows for the attractive computational properties of the spectral-element to be preserved, while allowing for the incorporation of attenuation in both fluid and solid regions.