Inconsistency of tissue-specific alpha values

[jkosciessa]

Hi all,
I noticed inconsistent comments in the config.yaml, and a potential source for errors when setting the attenuation coefficient. Despite presumed dependence on the fundamental frequency, values are currently hardcoded and relate to 250kHz fundamental frequency (cohering with source_freq_hz) instead of the specified default_freq of 500k. These parameters can be overwritten by a study-specific script, but this seems like a pretty significant source of potential deviations in acoustic intensity and heating. Would be best to always calculate frequency-specific values, or at least document this as a specification to look out for.

Here is @eleonoracarpino's original comments in our internal channels (12th Oct. 2023):

hello simulations people! we have noticed that the skull alpha coefficients that we have been using in the simulations are not a really good fit. therefore, from now on we would suggest to calculate your skull alpha coefficients as follow:

medium.alpha_power = 2 
medium.alpha_coeff = 13.3 * (f [MHz]) / (f [MHz])^2 

where f [MHz] is the fundamental frequency expressed in MHz, so either 0.5 or 0.25 for our transducers
some info: 13.3 dB/ (MHz . cm) is a good estimate of the attenuation across the whole skull bone (assuming y = 1) for all frequencies. and we choose y = 2 for the simplicity of the calculations in k-wave. This will effect the simulations giving us lower values of isppa in the brain, because for example for the 250 kHz transducer the new medium.alpha_coeff = 53.2, while before we had 4 or 8. notice also that this value is unique for the whole skull and does not allow us to distinguish between cortical and trabecular bone; we found out talking with Brad; they don't have a good estimate for cortical vs trabecular but they have a general estimate for the skull (which is easier to get because you just send a wave trough the skull and measure what it remains on the other side of the skull). they also run their simulations with a unique value (and they are in parallel developing the pseudoCT to overcome this). The assumption is that 13.3 is a good estimate for all frequencies. Therefore the relationship 13.3 * f^1 = alpha_coeff * f^2 holds; for the 500 kHz transducer, new medium.alpha_coeff = 26.2 according to this formula, but other methods are possible in this case because there is more data available

[jkosciessa]

Just noticed that the updated values cited above are also a feature of the updated branch that merged @eleonoracarpino's changes. Both setups (see updated default.yaml waiting to be merged) would be inconsistent with our target settings as described in the comment above.

[jkosciessa]

The current comment explains alpha_0_true as
alpha_0 in alpha = alpha_0*frequency^alpha_power.

Not sure if this is true, or whether this directly specifies the alpha_coeff .

Update: it is used in two functions:

• in get_alpha_coeff:
  alpha_db_cm = medium.alpha_0_true*((freq/1e6).^medium.alpha_power_true)
• in setup_medium: fit in fitPowerLawParamsMulti: a = a0*f^y

[jkosciessa]

As per the example above:

freq = 500000
medium.alpha_0_true = 13.3/(freq/1e6); % = 26.6 | refered to as alpha_coeff above; *(f [MHz]) / (f [MHz])^2 can be simplified into  /f [MHz]
medium.alpha_power_true = 2
alpha_db_cm = medium.alpha_0_true*((freq/1e6).^medium.alpha_power_true)

Thus, alpha_db_cm will scale linearly according to fundamental frequency.

Note that in the previous default_config, medium.alpha_power_true = 1.2 was specified in the comment for 250kHz transducers. I do not know where this value came from.

TLDR: medium.alpha_0_true=26.6 for 500kHz applications following the logic on top.

[jkosciessa]

Commit fd82fb4 adjusts the default config to the 500kHz parameters described above.

Instructions of what parameters to adjust for other fundamental frequencies (e.g., 250/300 kHz) remains largely missing.

[jkosciessa]

I wonder whether there is a significant bug regarding the alpha values. In multiple parts of the current code, the medium mask creation seems to not extract the proper alpha_0 value from the config, but rather alpha_power_true (e.g., see here). Given that across tissue types, alpha_power_true=1.2 in the past and based on the rationale above alpha_power_true=2 now, this would intuitively mean that absorption was overestimated in skin and brain, and underestimated in the skull.

[jkosciessa]

Here is the table with the ITRUSST benchmarks. Notably, these are the frequency-specific alphas at 500 kHz. At 1 MHz, values should be higher. Based on Eleonora's discussion above alpha should be updated for both cortical and trabecular bone.

Bone (overall): 13.3 dB/cm at 1 MHz [Pinton et al., 2012]
500 kHz: interpolated between 6.6 (13.3 x 0.5^1) and 3.3 (13.3 x 0.5^2)

0.3–1 Np/cm = 2.6 - 8 dB/cm [500kHz, Hüter et al., 1952]
0.7–1.2 Np/cm = 6 - 10 dB/cm [500kHz, Webb et al., 2021]

References

Pinton, Gianmarco, et al. "Attenuation, scattering, and absorption of ultrasound in the skull bone." Medical physics 39.1 (2012): 299-307.

T. F. Hüter, “Messung der Ultraschallabsorption im menschlichen
Schädelknochen und ihre Abhängigkeit von der Frequenz,” Die Naturwissenschaften, vol. 39, no. 1, pp. 21–22, 1952.

Webb, T. D., Leung, S. A., Ghanouni, P., Dahl, J. J., Pelc, N. J., & Pauly, K. B. (2021). Acoustic Attenuation: Multifrequency Measurement and Relationship to CT and MR Imaging. IEEE transactions on ultrasonics, ferroelectrics, and frequency
control, 68(5), 1532–1545. https://doi.org/10.1109/TUFFC.2020.3039
743

[jkosciessa]

Based on these values, the following alpha_0_true should be used in combination with alpha_power_true=1

f = 0.5 % 500 kHz [data reported in benchmarks]
alpha_power_true = 1
0.2/ f^alpha_power_true = 0.4 % skin
0.3/ f^alpha_power_true = 0.6 % brain
13.3 % skull [already 1 MHz estimate, approx. 4-8 at 500 kHz]

The skin and brain setup is largely consistent with the previous default config (except that they may not have been properly used to the bug noted above), with the only deviation being that alpha_power_true=1.2 in the past.

If the bug about reading in alpha_power_true instead of alpha_0_true (described above) is indeed as I expect, this would mean that skull attenuation may have been frequently underestimated by a factor of 10-20, while skin/brain attenuation may have been overestimated by a factor of 2.

Perhaps relatedly, the calculation of alpha_0_true for pseudoCTs (frequency-independent) currently yields alpha_0_true values of about 5-7, so about 4-5 times lower than the value that would be used for the binary skull segmentation based on the reasoning above. I am not sure where that equation comes from, but this should lead to unintuitive differences between setups. This relates to #43.

alpha_0_true(medium_masks==label_i) = 4+4.7.*[1-(pseudoCT_cropped(medium_masks==label_i)-1000 -HU_min)./(HU_max-HU_min)].^0.5;

[jkosciessa]

Based on Fig 2 in https://pubs.aip.org/asa/jasa/article/136/4/1499/905255/Modeling-power-law-absorption-and-dispersion-in, absorption coefficients should increase with fundamental frequency. My current results rather double the absorption for 500 kHz.

Edit: Nevermind, I became confused by the processing steps. The values above are the alpha_0 that is fitted when using a fixed alpha exponent of 2 in the function fitPowerLawParamsMulti; see it's call here. This should be correct:

53.2*0.5MHz^2 = 13.3

So I think my updated setup may work as expected now (for the binary skull implementation; the deviation for pseudoCTs as described above remains unresolved; perhaps a scaling factor is needed there as the config expects values specified for 0.5 MHz, not 1 MHz).

[jkosciessa]

From Carpino et al. (2024). Transcranial ultrasonic stimulation of the human amygdala to modulate threat learning. MSc thesis.

The absorption coefficients in the updated branch would now match this reported specification for the skull, while skin and brain absorption coefficients would be double what is reported here, given that they have been specified in the default config for 500 kHz instead of 1 MHz. Note for instance that the default config specified alpha0 for an exponent of 1.2 (see here). Alpha0 is then internally always rescaled according to an exponent of 2 as reported in the table. But I am not sure whether this rescaling had been accounted for in the reported values.

The error of reading in a0 = y seems to also have affected the simulations Eleonora ran based on her codebase here.

[jkosciessa]

Based on my reading of the ITRUSST benchmarks at 500 kHz, the 500 kHz values may look like follows:

brain

0.6 x 0.5^1 = 1.2 x 0.5^2 = 0.3 [ITRUSST 0.5 MHz]

skin

0.4 x 0.5^1 = 0.8 x 0.5^2 = 0.2 [ITRUSST 0.5 MHz]

So the brain/skin values in Eleonora's table maybe should have been doubled when moving to a coefficient of 2. (Then again, this should be the 250kHz setup in her case.) For 500 kHz, the bone coefficient should perhaps also be halved, as 13.3 is the attenuation (alpha) at 1 MHz. Per her equation above:

bone

13.3 x 1^1 = 13.3 x 1^2 = 13.3 [Pinton et al., 1 MHz]
13.3 x 0.5^1 = 26.6 x 0.5^2 = 6.65 [interpolated, 500 kHz]
13.3 x 0.25^1 = 53.2 x 0.25^2 = 3.325 [interpolated, 250 kHz]

---

This fits the expected values at source_freq_hz of 500 kHz as fit in fitPowerLawParamsMulti.m.

---

What remains weird to me is that intuitively alpha0 should be a frequency-independent constant, with the frequency-dependency of alpha attenuation arising from the frequency-exponent term. This is also apparent from the plot above, whose values at 250 kHz would not match up with the calculated ones, as alpha0 now also changed with frequency.

[jkosciessa]

RE: pseudoCT, the alpha created based on the HU conversion, should already be the attenuation at 500 kHz based on page 3 of the supplementary material by Yaakub et al.; in @eleonoracarpino's implementation those would be turned into the alpha0 prefactors instead, and thus would depend on the absorption coefficient exponent 𝑦.

Given

alpha_db_cm = medium.alpha_0_true*((freq/1e6).^medium.alpha_power_true)

we need to convert the pseudoCT-based attenuation alpha into prefactors and exponents, e.g., for y=2:

alpha_0 = alpha_pseudoCT/(source_freq_hz[in Mhz]^y)

This step seems to be missing from Yakuub et al.'s implementation as well however. Intuitively, this should underestimate the attenuation given that the attenuation refers to values measured at 500 kHz and thus would understimate attenuation at 1 MHz by an exponential y (i.e., underestimated by a factor larger than 2). They could only straightforwardly be used as a prefactor at source_freq_hz = 1 MHz however (denominator = 1).

I suppose it should always be

alpha_0 = alpha_pseudoCT/(0.5^y)

but then the exponent really matters when computing solutions for any other central frequency beyond 500 kHz. For 500 kHz, this should lead to good convergence with the binary skull alpha = 6.65 above as a midpoint between 4 and 8.7 based on pseudoCT values.

26.6 = 6.65/0.5^2

[jkosciessa]

Just a bit more background info on how the parameters work within the pipeline:

setup_medium will convert the specified alpha prefactors and exponents to align with a fixed exponent and the specified sound speed. The function fitPowerLawParamsMulti can provide an optional output of the effective tissue attenuation alpha.

In theory, as long as prefactors and exponents in the config are specified to result in the desired attenuation at the specified fundamental frequency, this conversion should not have too much of an impact. I am not yet convinced that it doesn't have an impact in practice, e.g., by affecting the stability of solutions.

[KTZ228]

Thank you for diving into this to such an extent! To ensure future users that are not familiar with these issues, however, it might be good to create separate blocks of medium properties (so medium.500KHz) within the default config for the different fundamental frequencies given that we don't expect most users to play around with this.
To limit the number of required changes you could just add a line in load_parameters that does something like parameters.medium = fprintf('parameters.medium.%.0fKHz', (parameters.transducer.source_freq_hz/1000)).
It would also be good to throw an error if the values have not yet been calculated for the given frequency (say a new member will use new transducers with an as-of-yet unused fundamental frequency)

[jkosciessa]

Sounds like a good feature improvement. Intuitively, it should also be possible to settle on one particular value that holds across frequencies. Ultimately, the formula already considers fundamental frequency (alpha = a x freq ^ b). The problem is that we then also have to have a valid-ish estimate of exponent b (which now is chosen arbitrarily to match the fundamental frequency of interest).

[KTZ228]

I wasn't aware that it might be possible to settle on a set of values that hold for all fundamental frequencies. In that scenario using one set for all frequencies would indeed be preferred. However, if we cannot find a way to not vary b for different frequencies the proposed solution could be a foolproof fit.

[jkosciessa]

@KTZ228 The following may serve as a general definition and to a large extent converges with the definition used within k-Plan:

medium:
  water:
    sound_speed: 1500 # [m/s] ITRUSST benchmarks 
    density: 994 # [kg/m^3] Tissue Properties DB or waterDensity(37) function in kWave
    alpha_0_true: 2.17e-3  # k-Plan https://dispatch.k-plan.io/static/docs/simulation-pipeline.html#evaluating-plans
    alpha_power_true: 2 # k-Plan https://dispatch.k-plan.io/static/docs/simulation-pipeline.html#evaluating-plans
    thermal_conductivity: 0.60 # [W/m/°C] Tissue Properties DB
    specific_heat_capacity: 4178 # [J/kg/°C] Tissue Properties DB
  skull: # cortical bone 
    sound_speed: 2800 # [m/s], ITRUSST benchmarks 
    density: 1850 # [kg/m^3], ITRUSST benchmarks
    alpha_0_true: 13.3 # Pinton et al., 2011 [note: overestimates heating]
    alpha_power_true: 1 # k-Plan https://dispatch.k-plan.io/static/docs/simulation-pipeline.html#evaluating-plans
    thermal_conductivity: 0.32 # [W/m/°C] Tissue Properties DB 
    specific_heat_capacity: 1313 # [J/kg/°C] Tissue Properties DB
  brain:
    sound_speed: 1546 # [m/s] Tissue Properties DB 
    density: 1046 # [kg/m^3], ITRUSST benchmarks or waterDensity(temp_0)
    alpha_0_true: 0.59 # k-Plan https://dispatch.k-plan.io/static/docs/simulation-pipeline.html#evaluating-plans
    alpha_power_true: 1.2 # k-Plan https://dispatch.k-plan.io/static/docs/simulation-pipeline.html#evaluating-plans
    thermal_conductivity: 0.51 # [W/m/°C] Tissue Properties DB 
    specific_heat_capacity: 3630 # [J/kg/°C] Tissue Properties DB
  skin:
    sound_speed: 1610 # [m/s], ITRUSST benchmarks 
    density: 1090 # [kg/m^3], ITRUSST benchmarks
    alpha_0_true: 0.4 # ITRUSST benchmarks
    alpha_power_true: 1
    thermal_conductivity: 0.37 # [W/m/°C] Tissue Properties DB 
    specific_heat_capacity: 3391 # [J/kg/°C] Tissue Properties DB
  skull_trabecular: # trabecular bone
    sound_speed: 2300 # [m/s], ITRUSST benchmarks 
    density: 1700 # [kg/m^3], ITRUSST benchmarks
    alpha_0_true: 13.3 # alt: 21.1 benchmark Aubry et al., 2022 [8 at 500 kHz]
    alpha_power_true: 1 # alt: 1.4 Robertson et al., PMB 2017
    thermal_conductivity: 0.32 # [W/m/°C] Tissue Properties DB 
    specific_heat_capacity: 1313 # [J/kg/°C] Tissue Properties DB
  skull_cortical: # cortical bone
    sound_speed: 2800 # [m/s], ITRUSST benchmarks 
    density: 1850 # [kg/m^3], ITRUSST benchmarks
    alpha_0_true: 13.3 # alt: 10.5 benchmark Aubry et al., 2022 [4 at 500 kHz]
    alpha_power_true: 1 # alt: 1.4 Robertson et al., PMB 2017
    thermal_conductivity: 0.32 # [W/m/°C] Tissue Properties DB 
    specific_heat_capacity: 1313 # [J/kg/°C] Tissue Properties DB

[jkosciessa]

The values used by k-plan are documented here.

[jkosciessa]

Frequency-specific calculation of attenuation/absorption

• The function get_alpha_coeff may not be in use (and may not have been used for some time).
• fitPowerLawParamsMulti is used in setup_medium to estimate alpha_coeff with a fixed exponent of 2. It considers the fundamental frequency. kspaceFirstOrder2D expects medium.alpha_coeff as [dB/(MHz^y cm)] (see doc), i.e., not frequency-specific.
• run_heating_simulations initially converts attenuation to the frequency-specific estimate.