A Python FDTD method algorithm for 1D planar acoustic wave propagation: Simulating high-frequency ultrasound in the brain and beyond

A Python FDTD method Algorithm for 1D Planar Acoustic Wave Propagation: Simulating High-Frequency Ultrasound in the Brain and Beyond Nuno A.T.C. Fernandes1[0000-0003-1848-726X]*, Ana Arieira1[0000-0003-4824-682X], Betina Hinckel2[0000-0003-3963-5241], Filipe Silva1[0000-0003-3596-3328], Ana Leal1[0000-0002-2844-491X], Óscar Carvalho1[0000-0002-9447-8739] 1 Center for microelectricomechanical systems (CMEMS), University of Minho, 4800-058 Guimarães, Portugal 2 Department of Orthopaedic Surgery, William Beaumont Hospital, Royal Oak, MI, USA [email protected] Abstract. Non-invasive techniques, such as high-frequency ultrasound, have emerged as promising therapeutic tools for neurological disorders, including Parkinson’s disease and Alzheimer’s disease. By targeting specific brain regions, ultrasound stimulation modulates neural activity and induces beneficial physiological responses. However, simulating high-frequency acoustic wave propagation in biological tissues presents computational challenges due to the high spatial and temporal resolution required to satisfy the low CourantFriedrichs-Lewy (CFL) condition for numerical stability and accuracy. This paper introduces a novel Python algorithm optimized for planar wave propagation, enabling efficient one-dimensional simulations of high-frequency ultrasound. Utilizing the finite difference time domain (FDTD) method, the algorithm incorporates material-specific properties, including density, sound speed, and frequency-dependent attenuation, to model heterogeneous tissue structures such as skin, bone, cerebrospinal fluid, and brain tissue. The method accurately captures key acoustic phenomena, such as impedance mismatching and wave reflection, facilitating detailed analysis of energy transmission and absorption in complex biological interfaces. The algorithm’s performance is compared with COMSOL Multiphysics, which is inherently limited to two and threedimensional acoustic wave propagation. By reducing the problem to one dimension, the proposed method simplifies computational complexity while preserving key wave interactions, enabling early-stage analysis with lower computational costs. Beyond biomedical applications, this approach is broadly applicable to any system governed by acoustic wave equations. By significantly reducing computational demands, it accelerates preliminary studies of wave propagation through multilayered media, contributing to the development of efficient ultrasound-based therapeutic models and advancing acoustic research. Keywords: Finite-Difference Time-Domain, Ultrasound Propagation, Biological Tissue Acoustics, Transcranial Ultrasound, Acoustic Attenuation, Neuromodulation, Simulation. 2 1 Introduction Ultrasound is a powerful modality for both medical diagnostics and non-destructive testing (NDT) due to its ability to penetrate layered media and provide real-time feedback on internal structures [1], [2]. Ultrasound is widely used in biomedical applications for imaging [3], therapy [4], and real-time monitoring of both soft and hard tissues [5], offering a non-invasive and versatile diagnostic and therapeutic tool. In industrial NDT, ultrasonic transducers are widely used in industrial non-destructive testing to detect internal flaws in metallic components, ensuring structural integrity without damaging the material [6], [7]. Despite the differences in application, both domains share common challenges related to acoustic wave propagation through heterogeneous, multilayered materials. Neurological disorders such as Alzheimer's disease, Parkinson's disease, and stroke are the leading cause of disability and the second leading cause of death globally, accounting for over 10 million deaths and 349 million disability-adjusted life years in 2019 [8], [9], [10]. Traditional diagnostic and therapeutic approaches for neurological disorders often rely on invasive methods, unraveling growing need for non-invasive, precisely targeted technologies such as transcranial magnetic stimulation and focused ultrasound, which offer safer alternatives for modulating neural activity and delivering treatments [11], [12]. Among non-invasive stimulation techniques, Ultrasoundbased techniques are emerging as promising non-invasive methods for brain stimulation [13], facilitating drug delivery across the blood-brain barrier [14], and enabling real-time neurological monitoring [15]. Accurately modeling the propagation of ultrasound through cranial structures, including the skin, skull, and underlying brain tissue is essential for optimizing the safety, precision, and efficacy of emerging neurological acoustic technologies. This research is driven by the broader objective of advancing non-invasive brain health diagnostics and therapeutic interventions through highfidelity simulation of acoustic wave behavior in anatomically and acoustically heterogeneous cerebral environments. Numerical modeling using finite-difference time-domain (FDTD) methods is widely recognized as a powerful approach for simulating acoustic wave propagation in inhomogeneous media [16], [17], [18]. FDTD enables direct computation of pressure and velocity fields with high spatial and temporal resolution and is particularly effective in incorporating spatially varying material properties, frequency-dependent attenuation, and complex boundary conditions. Studies have demonstrated the versatility of FDTD in modeling nonlinear effects, tissue-specific damping, and realistic source injection in both static and dynamic acoustic environments, reinforcing its relevance in biomedical and engineering acoustics applications [19], [20]. In this study, we present a one-dimensional FDTD model of ultrasound propagation through a simplified human head cross-section, including a piezoelectric transducer, a polymeric matching layer, and sequential layers of skin, bone, and brain tissue. The model incorporates frequency-dependent attenuation, impedance mismatches, and spatial filtering to emulate ideal transducer configurations. The goal is to provide a lightweight yet realistic simulation tool for investigating acoustic behavior in 3 biological conditions, which can support both medical device development and fundamental biophysics research. 2 Theoretical background In one dimension, the propagation of longitudinal acoustic waves in a medium is governed by two coupled first-order equations: the momentum equation derived by Newton’s first law (Eq. (1)) and the continuity equation derived by mass conservation (Eq. (2)) [21]. ∂c/∂t =-∂p/(ρ·∂x) (1) ∂p/∂t =-K·∂c/∂x (2) Together, these equations define a coupled system that supports wave propagation. Their linearity assumes small-amplitude oscillations and neglects nonlinear, thermal, or viscous effects, making them particularly suitable for modeling diagnostic or therapeutic ultrasound under moderate intensity conditions. When an acoustic wave encounters a change in material properties, part of the wave is reflected, and part is transmitted. The reflection and transmission are determined by the acoustic impedance, which can be simplified for planar waves to Eq. (3). Z =ρc (3) The mismatch of impedances can limit the amount of pressure that is transmitted. As such, the reflection and transmission of a planar acoustic wave can be described by the reflection and transmission coefficient presented in Eq. (4) and Eq. (5) respectively. R =(Z2-Z1)/(Z2+Z1) (4) T =2·Z2/(Z2+Z1) (5) Acoustic attenuation in soft tissues and bone is frequency-dependent and is typically characterized in units of dB/MHz/cm. To incorporate this into the simulation, attenuation coefficients Np/m are calculated as Eq (6). α=(αdB/(20·log10e))· (f/106)·100 (6) This attenuation arises from a combination of viscous, thermal, and structural relaxation processes within biological tissues, which convert a portion of the acoustic energy into heat [22]. The frequency dependence of this phenomenon is especially relevant in biomedical applications, where higher frequencies result in greater energy loss, directly affecting the amplitude, waveform shape, and penetration depth of the acoustic signal. Since these dissipative effects are not included in the theoretical derivation of the linear acoustic wave (Eq. (1) and Eq. (2)), attenuation is introduced as an 4 ad hoc correction during the numerical implementation within the FDTD update steps. 3 Methodology A first-order acoustic FDTD solver was implemented in Python, using explicit time integration and a staggered grid. Particle velocity (Eq. (7)) and pressure (Eq. (8)) are defined at alternating spatial positions to improve numerical stability and reduce dispersion. cin+1/2= cin-1/2-(pi+1n-pin)·Δt/(ρi·Δx) (7) pin+1/2= pin -(ci+1n+1/2-cin+1/2)·Δt·Ki/Δx (8) This formulation inherently accounts for reflection and transmission at material interfaces, as the spatial variation of density and bulk modulus directly governs wave behavior across boundaries. Material parameters were spatially assigned based on literature values [23], [24], [25], and the transducers present in the laboratory, as follows in Table 1. Table 1. Material parameters used for simulation. Material c (m/s) ρ (kg/m3) α (dB/MHz/cm) Skin 1540 1100 1.0 Bone 4080 1900 20.0 Brain 1560 1040 0.8 Transducer (PZT) 4200 7700 0 Matching Layer 2700 1500 0 The simulation domain spans 18 cm, consisting of a 6 cm brain-side buffer, a 3 cm biological interest region, and a 9 cm transducer-side buffer to avoid reflections. The transducer and matching layer were dimensioned according to their central frequency. The central 3 cm region includes skin, bone, and brain layers, each with user-defined thicknesses and properties, according to the literature [26]. Several regions of the brain were studied, and their respective thicknesses are presented in Fig. 1. The wave was generated using a single-cycle sine burst modulated by a Hann window to reduce spectral leakage. The waveform was injected into the particle velocity field to minimize numerical instability. The amplitude was scaled based on the transducer impedance to match a target pressure value. Five distinct frequencies were implemented: 0.5, 1, 2, 4 and 5 MHz. Biological tissue attenuation was implemented as an exponential decay applied directly to the pressure field during each update step. Attenuation coefficients were converted from dB/MHz/cm to Nepers/m and scaled by local speed of sound and time step. The attenuation is implemented as an exponential decay factor per time step (Eq, (9)), which damps the pressure field (Eq. (10)). 5 Fig. 1. Anatomical map showing skin and skull thicknesses across key cranial regions relevant for transcranial ultrasound targeting. γ =α·c·Δt (9) pin= pin·γ (10) To ensure numerical stability and accuracy, the spatial resolution Δx was chosen based on the shortest wavelength in the system Δx≤λmin/10. The time step was set according to the Courant–Friedrichs–Lewy condition Δt≤CFL·Δx/cmax, with CFL=0.1. Time-domain pressure signals were recorded at several positions: at the transducerskin barrier, skin-bone barrier, bone-brain barrier, and 5, 10 and 15 mm brain depth. To enhance signal visibility and suppress numerical noise, low-amplitude values (<107 Pa) were clipped, as they will be hard to detect with typical instrumentation. 4 Results The results of the simulations of the Python algorithm and COMSOL simulations are displayed in Fig. 2. Our findings clearly demonstrate the expected frequency-dependent attenuation of acoustic pressure as the wave propagates through brain tissues. Higher frequencies (4–5 MHz) exhibit a much steeper decline in acoustic pressure across all regions, while lower frequencies (500 kHz–1 MHz) maintain relatively higher amplitudes even at greater depths. This aligns with theoretical predictions and experimental data, as biological tissue attenuation increases approximately linearly with frequency. The comparison between the Python FDTD simulations and COMSOL outputs shows 6 strong agreement across all frequencies and anatomical regions, validating the numerical implementation and confirming that both attenuation and reflection/transmission effects due to acoustic impedance mismatches are accurately modelled. Fig. 2. Pressure attenuation across brain depth for six cranial regions and multiple frequencies, illustrating frequency-dependent losses and regional variability in acoustic transmission. It is also possible to observe the spatial variation in attenuation across different cranial regions, reflecting anatomical differences such as skull and skin thickness. For instance, the mastoid and temporal regions, which have thinner skull sections, exhibit less attenuation, whereas the frontal and vertex regions show steeper drops in pressure, consistent with denser bone and layered impedance. Since these values exclude the propagation path through the skin and skull, it's important to recognize that a significant portion of the acoustic signal is already attenuated before it even reaches brain tissue. The skin and skull layers, especially the cortical bone, present substantial acoustic impedance mismatches and frequency-dependent attenuation, which result in reflection, scattering, and absorption of the incoming wave. As a result, the pressures plotted at 0 mm brain depth do not represent the original transducer output but rather a diminished signal that has already undergone considerable energy loss. 7 Ultrasound neuromodulation holds significant therapeutic potential for a range of neurological conditions, particularly when its ability to penetrate different cranial regions is matched with the anatomical region of the disease. Alzheimer’s disease, a progressive neurodegenerative disorder marked by memory loss and cognitive decline, primarily affects the temporal and frontal lobes [27], these deeper regions benefit from lower-frequency ultrasound (e.g., 500 kHz–1 MHz), which offers better penetration. Frontotemporal dementia, which alters behaviour, personality, and language, also targets these frontal and temporal structures [28] and may similarly benefit from low-frequency neuromodulation. Posterior cortical atrophy, a visual variant of Alzheimer’s disease, affects the occipital lobe, which is relatively superficial [29] and thus accessible to higher-frequency ultrasound that offers finer spatial resolution. Similarly, multiple sclerosis, a chronic autoimmune disease that damages the protective covering of nerve fibers (myelin), often affects multiple, diffusely located brain areas including the vertex region [30], demanding careful targeting strategies. Finally, schizophrenia and major depressive disorder, which involve cognitive dysfunction and emotional dysregulation, are strongly associated with the prefrontal cortex, a relatively superficial region, making it an ideal target for transcranial ultrasound interventions [31]. Adapting ultrasound parameters to match both the disease characteristics and the anatomical constraints of each region reinforces the promise of neuromodulation as a precise, non-invasive therapeutic tool. Fig. 3. Comparison of simulation runtimes for Python and COMSOL across frequencies. The comparison of simulation runtimes between Python and COMSOL illustrated in Fig. 3 clearly demonstrates the computational efficiency of the Python-based FDTD implementation, particularly as frequency increases. While both methods perform relatively quickly at lower frequencies, the runtime for COMSOL grows exponentially with frequency, reaching over 3.57 hours at 5 MHz, due to the increased spatial resolution and mesh refinement required for higher-frequency wave propagation in finite element models. In contrast, the Python solver maintains a much smaller computational footprint, scaling more gracefully with frequency thanks to its explicit time-stepping and structured grid approach. This efficiency makes the Python imple- 8 mentation highly suitable for iterative studies, parametric sweeps, or real-time applications, especially when rapid feedback is essential. The trade-off, of course, is that COMSOL offers greater modelling flexibility and precision for complex geometries, but for structured 1D or layered media problems like transcranial wave modelling, Python provides an excellent balance of speed and fidelity. 5 Conclusion This study presented a computational framework for modelling ultrasound wave propagation through layered cranial tissues using a 1D FDTD method. By incorporating spatially varying acoustic properties, frequency-dependent attenuation, and anatomical variability in skin and skull thickness, the model effectively captures the complex interactions that influence transcranial acoustic transmission. Validation against COMSOL simulations demonstrated strong agreement in pressure profiles across multiple cranial regions and frequencies, while also highlighting the substantial computational efficiency of the Python implementation. The results confirm the critical role of frequency selection in balancing spatial resolution and penetration depth, with lower frequencies enabling deeper brain access and higher frequencies offering more localized energy delivery. Regional differences in attenuation further emphasize the importance of personalized acoustic targeting, particularly in the context of neuromodulation therapies for conditions such as Alzheimer’s disease, epilepsy, and depression. Overall, this lightweight simulation approach offers a valuable tool for exploring and optimizing non-invasive ultrasound applications in brain health monitoring and intervention. 6 Acknowledgements This work was funded by the National Foundation for Science and Technology of Portugal (FCT) under the project “BrainStimMap – Mapping and modelling the transmission profile of optomechanical waves in the brain to optimize transcranial stimulation against brain disorders" with the reference PTDC/EME-EME/1681/2021 and the individual PhD grant with reference 2022.11063.BD. References [1] C. M. I. Quarato et al., ‘A Review on Biological Effects of Ultrasounds: Key Messages for Clinicians’, Diagnostics, vol. 13, no. 5, p. 855, Feb. 2023, doi: 10.3390/diagnostics13050855. [2] A. Akundi, T.-L. Tseng, M. F. Rahman, and E. Smith, ‘Non-Destructive Testing (NDT) and Evaluation Using Ultrasonic Testing Equipment to Enhance Workforce Skillset for Modern Manufacturing’, in 2018 ASEE Annual Conference & Exposition Proceedings, Salt Lake City, Utah: ASEE Conferences, 2018. doi: 10.18260/1-2--30842. 9 [3] G. Cloutier, F. Destrempes, F. Yu, and A. Tang, ‘Quantitative ultrasound imaging of soft biological tissues: a primer for radiologists and medical physicists’, Insights Imaging, vol. 12, no. 1, Dec. 2021, doi: 10.1186/s13244-021-01071-w. [4] D. L. Miller et al., ‘Overview of Therapeutic Ultrasound Applications and Safety Considerations’, Journal of Ultrasound in Medicine, vol. 31, no. 4, pp. 623–634, Apr. 2012, doi: 10.7863/jum.2012.31.4.623. [5] C. Wang et al., ‘Bioadhesive ultrasound for long-term continuous imaging of diverse organs’, Science, vol. 377, no. 6605, pp. 517–523, Jul. 2022, doi: 10.1126/science.abo2542. [6] ‘Review of Ultrasonic Testing for Metallic Additively Manufactured Parts’, in Additive Manufacturing Design and Applications, ASM International, 2023, pp. 310–323. doi: 10.31399/asm.hb.v24a.a0006982. [7] A. Chabot, N. Laroche, E. Carcreff, M. Rauch, and J.-Y. Hascoët, ‘Towards defect monitoring for metallic additive manufacturing components using phased array ultrasonic testing’, J Intell Manuf, vol. 31, no. 5, pp. 1191–1201, Jun. 2020, doi: 10.1007/s10845-019-01505-9. [8] V. L. Feigin et al., ‘The global burden of neurological disorders: translating evidence into policy’, The Lancet Neurology, vol. 19, no. 3, pp. 255–265, Mar. 2020, doi: 10.1016/s1474-4422(19)30411-9. [9] V. L. Feigin et al., ‘Global, regional, and national burden of neurological disorders, 1990–2016: a systematic analysis for the Global Burden of Disease Study 2016’, The Lancet Neurology, vol. 18, no. 5, pp. 459–480, May 2019, doi: 10.1016/s1474-4422(18)30499-x. [10] C. Ding et al., ‘Global, regional, and national burden and attributable risk factors of neurological disorders: The Global Burden of Disease study 1990–2019’, Front. Public Health, vol. 10, Nov. 2022, doi: 10.3389/fpubh.2022.952161. [11] G. Chen et al., ‘Non-invasive brain stimulation effectively improves post-stroke sensory impairment: a systematic review and meta-analysis’, J Neural Transm, vol. 130, no. 10, pp. 1219–1230, Oct. 2023, doi: 10.1007/s00702-023-02674-x. [12] M. Bange, R. C. G. Helmich, A. A. Wagle Shukla, G. Deuschl, and M. Muthuraman, ‘Non-invasive brain stimulation to modulate neural activity in Parkinson’s disease’, npj Parkinsons Dis., vol. 11, no. 1, Apr. 2025, doi: 10.1038/s41531-025-00908-1. [13] G. Darmani et al., ‘Individualized non-invasive deep brain stimulation of the basal ganglia using transcranial ultrasound stimulation’, Nat Commun, vol. 16, no. 1, Mar. 2025, doi: 10.1038/s41467025-57883-7. [14] M. Gionso et al., ‘Ultrasound guided blood brain barrier opening using a diagnostic probe in a whole brain model’, Sci Rep, vol. 15, no. 1, Mar. 2025, doi: 10.1038/s41598-025-94660-4. [15] K. R. Murphy et al., ‘Optimized ultrasound neuromodulation for non-invasive control of behavior and physiology’, Neuron, vol. 112, no. 19, pp. 3252-3266.e5, Oct. 2024, doi: 10.1016/j.neuron.2024.07.002. [16] G. Petris, M. Cianferra, and V. Armenio, ‘A numerical method for the solution of the threedimensional acoustic wave equation in a marine environment considering complex sources’, Ocean Engineering, vol. 256, p. 111459, Jul. 2022, doi: 10.1016/j.oceaneng.2022.111459. [17] G. V. Norton and J. C. Novarini, ‘Finite-difference time-domain simulation of acoustic propagation in dispersive medium: An application to bubble clouds in the ocean’, Computer Physics Communications, vol. 174, no. 12, pp. 961–965, Jun. 2006, doi: 10.1016/j.cpc.2006.01.003. [18] M. Saeki et al., ‘FDTD simulation study of ultrasonic wave propagation in human radius model generated from 3D HR-pQCT images’, Physics in Medicine, vol. 10, p. 100029, Dec. 2020, doi: 10.1016/j.phmed.2020.100029.