Designing sensor input, cluster code, etc..

As a quick overview driving the following questions, we stimulated a certain brain area using a 0.5 Mhz single element focused transducer. The stimulation was pulsed, lasting for 500 ms. Another group simulated this setup using k-wave in this paper: Numerical evaluation of the skull for human neuromodulation with transcranial focused ultrasound
Jerel K Mueller, Leo Ai, Priya Bansal and Wynn Legon.

The basis of our questions stems from trying to essentially replicate the simulations they did in K-wave. We have thus far written a substantial amount of code trying to do it, but still have several questions from the experts. So, first, thanks for taking the time to help.

I am currently trying to simulate the pressure field produced by an ultrasound transducer when the ultrasound signal is passed through a human skull immersed in water. The transducer has the following specifications:
Focal length = 36 mm Diameter of aperture of transducer = 30 mm Radius of curvature of transducer = 30 mm Center frequency = 0.5 MHz

The skull is immersed in water and the transducer is placed above the skull with a thin layer of water between the skull and transducer i.e. the transducer isn’t placed directly on the skull. I am trying to record the pressure field at and around the focus which is on the other side of the skull as compared to the transducer. The transducer is driven by an input signal with the following specifications:
Acoustic frequency = 0.5 MHz Pulse duration = 360 μs Pulse repetition frequency = 1.0 kHz
Stimulus duration = 0.5 s Number of pulses = 500

I would be happy to post our code if that helps, or if you have code available (for example, from your paper: "Accurate simulation of transcranial ultrasound propogation for ultrasonic neuromodulation and stimulation). That paper of yours and the Legon one mentioned above described exactly what we are trying to do simulate.

1. The first problem we encounter is that MATLAB usually crashes while executing the code, which seems to be a memory issue. We have access to a HPC and a gpu-cluster. Do you have code available to run the rather large simulations on a GPU cluster? I believe the cluster is 8 Tesla V100 gpus. We also have access to more traditional multi-node clusters with 128+ cores. Alternatively, is there a simple way to split up the calculations on a local computer using parallel processing and recombine the sensor data afterwards?

2. A rather naive question (due to lack of electronics knowledge): Question is about replicating our physical setup.. the transducer was connected to a 40 watt linear amp with a 50dB gain, which was excited with a pulsed 2volt signal. Would this be properly implemented by changing the input signal to the transducer simply as an amplitude gain? Is there a straightforward way to set this up using the input? This is identical to the setup used in this paper: Transcranial focused ultrasound modulates the activity of primary somatosensory cortex in humans: Legon et al., Nature Neuroscience (2014).

3. Since we are trying to simulate a single-element spherically focused transducer having a focal length of 30 mm, would the radius of curvature also be 30 mm in the code?

4. Finally, is it better to set the source using source.p0 or source.p_mask for these types of simulation? Since what we are doing is a close approximation to the "Simulating Transducer Field Patterns Example", I assume using source.p_mask and source.p is better. Then, I believe this relates to my question 2 in how to make sure these values are being properly set.

(1) Yes, we have both an MPI-code (stable now for several years) and a multi-GPU code (still somewhat experimental). We are planning on releasing the MPI code later this year. Does the cluster you work with have singularity installed? Alternatively, do you (or someone you work with) have experience compiling libraries and binaries for linux?

(2) I think it's difficult to map directly from the electrical signal to an acoustic output. If possible, I'd suggest performing a characterisation measurement. There are at least two ways to go about this:

(a) Experimentally measure a 2D plane in water, and then derive a source plane which you input to k-Wave. This is explained in this paper. The codes to do this will be in the next release of k-Wave, but I'm happy to email them to you in advance.

(b) Experimentally measure the pressure at the focus in water. You can then use the formula for the focusing gain of a focused bowl transducer, and work backwards to calculate what pressure you need to set in k-Wave. Let me know if you need references for this.

(3) Yes.

(4) Use a source.p_mask (set to the shape of your transducer, e.g., using makeBowl) and source.p set to a time varying sinusoid at 0.5 MHz (e.g., using createCWSignals). To make the calculations easier, I would simply set the simulation time in k-Wave to be long enough to reach steady state (i.e., ignore the actual stimulus duration).

For 1.) the cluster does have singularity, from what I know. We have people whom have experience compiling libraries etc. I’m guessing the gpu cluster variant will be the fastest code for 3D simulations, correct? If you would send us the code, and any tips on making either the mpi or cuda work would be greatly appreciated.

For 2.). We are doing a water tank test right now to characterize the transducer and free water and thru a skull cap.
The code for the source plane (option a) would be great as well to get. For option b, yes, can you forward the references for backwards calculating the pressure that is needed.

My email is jmfine@asu.esu. But a google drive or Dropbox for those works well, too. Thank you again for all the help.

I am a graduate student working with Dr. Justin Fine for carrying out acoustic simulations involving the human skull as described in the previous posts. Please could you help me with the following:

1. Please could you email me the code for deriving a source plane from measured 2D plane in water for an acoustic transducer and the references of the formula for focusing gain of a focused bowl transducer to know what pressure to set in the K-wave simulation.
2. How do I approximate the time required for the simulation to reached steady state?

Regarding steady state, it depends on the medium heterogeneity. I would suggest running at least one simulation for a long time (e.g., 3 round trips of the domain length), and recording the pressure time series somewhere near the region of interest. You can then tell pretty easily when the transient behaviour has died out.