A method for optimizing magnetic resonance guided transcranial ultrasound neuromodulation target focusing

Abstract

The present application relates to the technical field of energy-based adaptive focusing, and discloses a method for optimizing magnetic resonance guided transcranial ultrasound neuromodulation target focusing. In combination with the iterative technique of the UTE image and the phase inversion algorithm, the short T2 water signal information in the skull captured by the UTE image is used to provide prior knowledge to reduce the number of calibration sequences; on the other hand, the DnCNN technology is also selected to resist the influence of noise, and the two can be well combined on the prDeep algorithm, so that the precise focusing of the transcranial ultrasound is finally realized.

Description

Technical Field

[0001] This invention relates to the field of energy-based adaptive focusing technology, specifically a method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing. Background Technology

[0002] In the fields of acoustic or optical imaging systems, the propagation of waves through scattering media has always been a fundamental physical problem. In acoustics, focusing ultrasound through the skull onto deep brain regions has become an emerging neuromodulation technique in recent years. It can non-invasively reach deep brain regions associated with neurological disorders (essential tremor, Parkinson's tremor, etc.) and focus the ultrasound within a few centimeters. However, the significant differences in sound velocity and strong attenuation between the skull, brain tissue, and degassed water lead to substantial distortion of the sound field at the target location, making skull-induced phase distortion a significant challenge.

[0003] Currently, numerical simulation is the only widely used clinical approach. It utilizes CT scans to acquire information about the patient's skull and builds a three-dimensional acoustic model. It mitigates these challenges by estimating the pressure amplitude at the focal point, the risk of skull heating, the size and location of the focal region, and heating within that region. However, the accuracy of three-dimensional simulation is limited by precision and relative positioning (the coordinates of the probe and patient's head in MR space must be repositioned for each treatment), and the full 3D acoustic simulation takes a long time to compute, making it unsuitable for real-time or rapid-response treatments, such as those for stroke.

[0004] In recent years, thanks to the continuous upgrades and development of MRI (Magnetic Resonance Imaging) hardware and ARFI (Acoustic Radiation Force Impulse) sequences, energy-based adaptive focusing technology has returned to the public eye. This technology uses reference acoustic intensities provided by MR-ARFI images to iteratively calculate the transmission matrix between the source and the target region, thereby using phased array ultrasound to compensate for phase loss caused by the skull, achieving transcranial focusing. For example, adaptive focusing methods based on Hadmaed, Zernike, and Random bases obtain corresponding displacement maps by inputting different encoding modes, thus estimating the phase difference caused by the transducer acoustic field and the ultrasound propagation medium to address the phase difference challenge caused by the skull. This technology does not require prior information from CT scans and is not limited by relative position. However, since iterative measurement of the transmission matrix requires the participation of ARFI sequences, the number of calibration sequences and the influence of image measurement noise are key factors in whether this technology can quickly and accurately measure the human skull transmission matrix, but current technologies cannot effectively address both simultaneously. Summary of the Invention

[0005] To address the aforementioned problems, the present invention aims to provide an optimized method for focusing transcranial ultrasound neuromodulation targets under magnetic resonance imaging (MRI) guidance. This method simultaneously solves the problems of excessive calibration cycles and experimental measurement noise; it also eliminates the need to scan the transducer and skull images before each treatment, overcomes the drawbacks of excessive iterations and susceptibility to noise, and achieves rapid and accurate transcranial focusing. The technical solution is as follows: A method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing includes the following steps: Step 1: Collect subject information: Scan the subjects with UTE and T1 sequences on MRI to obtain dicom images of the skull and brain, respectively; Step 2: Ultrasound simulation preparation: Input the skull and brain DICOM images into the K-Wave toolkit of acoustic simulation software in Matlab, use the regression curves of UTE and CT to calculate the skull porosity, and then derive the skull acoustic parameters, establish an acoustic model, and delineate the treatment area. Step 3: Preliminary simulation calculation: Perform time-reversal simulation on the focused treatment area on the computer. Calculate the input matrix by propagating ultrasound backwards from the target point. This is called a time-reversal mirror, or time delay; Step 4: Determine the focus position: Input the time delay obtained in the previous step into the actual treatment scenario, and use the ARFI sequence to trigger the ultrasound. Observe the actual focusing situation. If the ultrasound is accurately focused on the treatment area, the formal treatment will begin; otherwise, proceed to the next step. Step 5: Measure the transmission matrix: Using the acoustic model established in Step 2, calculate the transmission matrix value at the focal point. ; Step 6: Randomly input and collect transmission matrix information: Set the input matrix E to a size of 3. N × N A random matrix, where N is the number of ultrasonic array elements; input matrix Lines 1 to 10 The acoustic intensity values ​​at the target point after these 10 ultrasonic calibrations were recorded using ARFI sequences. ; Step 7: Iteratively calculate the new transfer matrix: , and Substituting the values ​​into the prDeep algorithm, we can calculate the conjugate of the new transmission matrix estimate, i.e., the time-reversed solution. Step 8: Verify if the new transfer matrix is ​​close to the true value: Input the obtained time-reversal solution into the actual transducer and use ARFI to observe whether the focal point is corrected. If so, the loop ends and the treatment begins; otherwise, return to step 6 and input... Lines 11-20 continue with ultrasonic calibration and ARFI recording of the target point acoustic intensity. Then, the new transmission matrix conjugate value is estimated by prDeep, and so on, until the focus is corrected.

[0006] Furthermore, the calculation of skull porosity in step 2 specifically involves: The skull is approximated as being made of water and bones of varying densities. An approximation of the porosity of the skull: ; In this case, the skull porosity map is directly linked to the Hausfield map: ; The acoustic properties of the skull include the speed of sound ,density and attenuation coefficient ,Right now: ; ; ; in, The linear attenuation coefficient of X-ray material at position x is the X-ray material at the tissue location. is the linear decay coefficient of water; The linear decay coefficient of the bone; This is the Heinz value for the corresponding position; The reference speed of sound for water; The reference sound velocity for the skeleton; The density of water; The density of bones; The minimum sound attenuation coefficient for the skeleton; It is the largest bone sound attenuation coefficient.

[0007] Furthermore, in step 5, the transmission matrix value at the focal point is calculated. H Specifically: Let N be the number of ultrasonic array elements, L be the number of points at different spatial locations in the region of interest, and f be the center frequency; when the nth vibration source emits a short ultrasonic pulse At time, at a point in space Some observed values ​​at [location] ; In the time domain, the relationship between the two is as follows: ; in, Represents convolution. For the values ​​of the transmission matrix elements; In the steady-state frequency domain, it is as follows: ; in, , and These are the Fourier transforms of the observed value P(t), the transmitted value H(t), and the input value E(t), respectively. yes A ×1 vector, representing The sound pressure amplitude of each unit within the target area during the second calibration; It is the first The frequency domain representation of the input matrix of the transducer array in the secondary calibration is ×1 vector; Represents the maximum number of measurements; yes A matrix, the frequency domain representation of the transmission matrix; Taking the conjugate transpose of the above equation, we get: ; Among them, superscript Represents the conjugate transpose, observation matrix Input matrix ; Given the conjugate of the input matrix ,Depend on Estimate each column The complex values ​​of each column are used to obtain the transmission matrix value at the focal point. .

[0008] Compared with the prior art, the beneficial effects of the present invention are: the present invention uses UTE images to provide prior information for the prDeep algorithm, thereby solving the problems of excessive calibration times and experimental measurement noise at the same time; and its advantage is that it does not need to scan transducer and skull images before each treatment to ensure that the relative positions remain unchanged during simulation, as is the case with numerical simulation methods. On the other hand, it also solves the shortcomings of previous iterative algorithms, such as excessive iteration times and large influence from noise. Finally, it can also achieve fast and accurate transcranial focusing effect. Attached Figure Description

[0009] Figure 1 The treatment flowchart is as follows: Dashed box 1 represents the operations performed under MRI during the actual treatment process; Dashed box 2 represents the cyclical process of adjusting the focus, where when the actual focus position is not in the treatment area and shifts using the ARFI sequence of MRI, the prDeep algorithm iteratively calculates the new transfer matrix value, thereby calculating a new input matrix to achieve iterative focusing until the actual focus position reaches the treatment area, at which point the cycle ends; Dashed box 3 represents the operations to be performed in the virtual computer environment, where MRI-UTE stands for Ultrashort Echo Time UTE Magnetic Resonance Imaging MRI; k-Wave is a MATLAB toolbox for time-domain simulation of acoustic fields; and prDeep is an algorithm for solving phase inversion problems. Detailed Implementation

[0010] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0011] To optimize the number of calibration sequences and reduce measurement noise, this invention proposes an iterative technique combining UTE (Ultra-short Time of Echo) images and a phase inversion algorithm. This technique utilizes the short T2 water signal information in the skull captured by the UTE images to provide prior knowledge and reduce the number of calibration sequences. Furthermore, DnCNN technology is added to combat noise. These two techniques are well integrated into the prDeep algorithm, ultimately achieving precise focusing of transcranial ultrasound. The treatment process is as follows: Figure 1 As shown.

[0012] The specific steps of the method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing in this invention are as follows: Step 1: Patients (subjects) are screened by clinicians and enter the treatment group. UTE and T1 sequences are scanned on MRI to obtain dicom (Digital Imaging and Communications in Medicine) images of the skull and brain, respectively.

[0013] Step 2: Input these images into the K-Wave toolkit of acoustic simulation software in Matlab, use the regression curves of UTE and CT to calculate the skull porosity, and then derive the skull acoustic parameters, establish an acoustic model, and delineate the treatment area.

[0014] Porosity calculations, based on a 2003 study by J.-F. Aubry et al., suggest that the skull can be approximated as being composed of water and bone of varying densities. Therefore, An approximation of the porosity of the skull: ; In this case, the skull porosity map is directly linked to the Hausfield map: ; To establish an acoustic model of the skull, it is also necessary to understand the acoustic properties of the skull, including density, sound velocity, and attenuation coefficient, i.e.: ; ; ; in, The linear attenuation coefficient of X-ray material at position x is the X-ray material at the tissue location. is the linear decay coefficient of water; The linear decay coefficient of the bone; Porosity of the skull; This is the Heinz value for the corresponding position; The reference speed of sound for water; The reference sound velocity for the skeleton; The density of water; The density of bones; The minimum sound attenuation coefficient for the skeleton; The maximum bone sound attenuation coefficient; It is the speed of sound; It is density; It is the attenuation coefficient.

[0015] Step 3: Perform a time-reversal simulation of the focused area on a computer to calculate the input at the transducer array, i.e., the time delay. Since the second-order partial derivatives of sound waves are time-invariant during propagation in the medium, we can calculate the input matrix E by propagating ultrasound backwards from the target point; this is called the time-reversal mirror (TRM).

[0016] Step 4: Input the time delay obtained in the previous step into the actual treatment scenario, and use ARFI sequence in conjunction with ultrasound to observe the actual focusing situation. Finally, determine whether it is focused on the target area. If it is, start the formal treatment (select different stimulation modes); otherwise, skip to the next step.

[0017] In actual treatment, the discrepancy between the ultrasound focal position and the acoustic simulation prediction results is due to the difference between the two transmission matrices. That is, the difference between the transmission matrix corresponding to the acoustic model and the transmission matrix of the subject's skull in the actual environment is too large, which leads to inaccurate focal prediction, i.e., focal shift.

[0018] Step 5: Using the acoustic model established in Step 2, calculate the transmission matrix value at the focal point. .

[0019] Step 6: Set the input matrix E to a size of 1. A random matrix (Bernoulli distribution +1, -1); input E Lines 1 to 10 The acoustic intensity values ​​at the target point after these 10 ultrasonic calibrations were recorded using ARFI sequences. .

[0020] Step 7: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation would require the full context.] , and Substituting the values ​​into the prDeep algorithm, we can calculate the conjugate of the new transfer matrix estimate, which is the time-reversed solution.

[0021] Step 8: Input the newly obtained time-reversal solution into the actual transducer. Use ARFI to observe whether the focal point has been corrected. If so, the loop ends and the treatment begins; otherwise, return to Step 6 and input... Lines 11-20 continue with ultrasonic calibration and ARFI recording of the target point acoustic intensity. Then, the new transmission matrix conjugate value is estimated by prDeep, and so on, until the focus is corrected.

[0022] In practice, the transmission matrix is ​​always unknown and should be determined using ultrasonic calibration by observing the ultrasonic output sound field through different input encoding modes. Therefore, the forward problem of measuring the transmission matrix can be standardized as: assuming the input ultrasonic waves are completely known, which transmission matrix best explains the currently observed output sound field?

[0023] Formally, let the number of ultrasonic array elements be N (N = 1:N), the number of points at different spatial locations in the region of interest (ROI) L be L (L = 1:L), and the center frequency be f = 250 kHz. When the nth vibration source emits a short ultrasonic pulse... At time, at a point in space The partial observations at that location (the square root of the local sound intensity) are: Because the transformation is linear, the relationship between the two in the time domain is: ,in, Represents convolution.

[0024] In the steady-state case (matrix form), it can be concisely written in the frequency domain, as shown below: ; The above equation can also be understood as the response of a monochromatic linear system, where, , and Observed values Transmitted values and input values Fourier transform; yes The ×1 vector represents the sound pressure amplitude of each unit within the target area; It is the input matrix of the transducer array in the μ-th calibration, and it is... ×1 vector; Represents the maximum number of measurements; yes The matrix is ​​called the transfer matrix; each element represents an input unit. With the target output point The transmission matrix contains not only information about the patient's skull, but also information about the relative positions of the transducer and the skull, as well as information about the free sound field. Therefore, in essence, the process of sound waves scattering through the skull, regardless of the complexity of the medium, can be fully described by the transmission matrix (TM).

[0025] Taking the conjugate transpose of the above equation, we get:

[0026] Among them, superscript Represents the conjugate transpose, observation matrix Input matrix .

[0027] The above equation is a "classic" phase inversion (pr) problem, that is, given the conjugate of the input matrix... ,So Each column is used for estimation The complex values ​​in each column.

[0028] As can be seen from the derivation of the above formula, the transcranial focusing task can be transformed into a phase inversion problem. The phase inversion problem arises when one wants to reconstruct a complex vector in a linear system from only the magnitude of a given measurement. However, traditional phase retrieval algorithms are difficult to implement in the presence of noise. Therefore, this invention uses the prDeep algorithm. This choice is based on two main advantages: 1. prDeep is created by utilizing a denoising regularization framework and a convolutional neural network denoiser, and has good noise resistance performance; 2. This general framework enables its use for calibrating TM, and the number of calibration sequences is related to the initial input complex value.

[0029] prDeep comprises two parts: prRED and DnCNN, exhibiting excellent adaptability and noise resistance. RED can be combined with any denoiser to regularize the general inverse problem. The process of using RED to solve the pr problem is called prRED. The other part involves the selection of a denoiser; DnCNN is a neural network used to remove Gaussian white noise from natural images. In this application, the framework we used includes four DnCNN networks with different noise levels, trained with 300,000 overlapping patches extracted from 400 images, adding Gaussian white noise, using the mean squared error between the noise-free actual image and the denoised reconstructed image as the cost function, and a final training rate of 0.00001.

Claims

1. A method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing, characterized in that, Includes the following steps: Step 1: Collect subject information: Scan the UTE and... on MRI. T 1 sequence, to obtain DICOM images of the skull and brain respectively; Step 2: Ultrasound simulation preparation: Input the skull and brain DICOM images into the K-Wave toolkit of acoustic simulation software in Matlab, use the regression curves of UTE and CT to calculate the skull porosity, and then derive the skull acoustic parameters, establish an acoustic model, and delineate the treatment area. Step 3: Preliminary simulation calculation: Perform time reversal simulation on the focused treatment area on the computer. By propagating ultrasound backwards from the target point, calculate the input matrix E, which is called the time reversal mirror, i.e., time delay; Step 4: Determine the focus position: Input the time delay obtained in the previous step into the actual treatment scenario, and use ARFI sequence in conjunction with ultrasound to observe the actual focusing situation. Finally, determine whether it is focused on the target area. If it is, start the formal treatment; otherwise, skip to the next step. Step 5: Measure the transmission matrix: Using the acoustic model established in Step 2, calculate the transmission matrix value at the focal point. H ; Step 6: Randomly input and collect transmission matrix information: Set the input matrix E to a size of 3. N × N random matrix, N The number of ultrasonic array elements; input matrix E Lines 1 to 10 E 10 The acoustic intensity values ​​at the target point after these 10 ultrasonic calibrations were recorded using ARFI sequences. P 10 ; Step 7: Iteratively calculate the new transfer matrix: P 10 , E 10 and H Substituting the values ​​into the prDeep algorithm, we can calculate the conjugate of the new transmission matrix estimate, i.e., the time-reversed solution. Step 8: Verify if the new transfer matrix is ​​close to the true value: Input the obtained time-reversal solution into the actual transducer and use ARFI to observe whether the focal point is corrected. If so, the loop ends and the treatment begins; otherwise, return to step 6 and input... E Lines 11-20 continue with ultrasonic calibration and ARFI recording of the target point acoustic intensity. Then, the new transmission matrix conjugate value is estimated by prDeep, and so on, until the focus is corrected.

2. The method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing according to claim 1, characterized in that, The calculation of skull porosity in step 2 is specifically as follows: The skull is approximated as being made of water and bones of varying densities. An approximation of the porosity of the skull: ； In this case, the skull porosity map is directly linked to the Hausfield map: ； Acoustic properties of the skull, including sound velocity ,density and attenuation coefficient ,Right now: ； ； ； in, The linear attenuation coefficient of X-rays at the corresponding tissue location; is the linear decay coefficient of water; The linear decay coefficient of the bone; This is the Henness value at the corresponding position; The reference speed of sound for water; The reference sound velocity for the skeleton; The density of water; The density of bones; The minimum sound attenuation coefficient for the skeleton; It is the largest bone sound attenuation coefficient.

3. The method for optimizing magnetic resonance-guided transcranial ultrasound neuromodulation target focusing according to claim 1, characterized in that, In step 5, the transmission matrix value at the focal point is calculated. Specifically: Assume the number of ultrasonic array elements N The number of points at different spatial locations within a region of interest in space. L and center frequency f ; When the nth vibration source emits a short ultrasonic pulse At time, at a point in space Some observed values ​​at [location] ; In the time domain, the relationship between the two is as follows: ； in, Represents convolution. For the values ​​of the transmission matrix elements; In the steady-state case, in the frequency domain, it is as follows: ； in, , and Observed values Transmitted values and input values Fourier transform; yes A ×1 vector, representing The sound pressure amplitude of each unit within the target area during the second calibration; It is the first The frequency domain representation of the input matrix of the transducer array in the secondary calibration is ×1 vector; Represents the maximum number of measurements; yes A matrix, the frequency domain representation of the transmission matrix; Taking the conjugate transpose of the above equation, we get: ； Among them, superscript Represents the conjugate transpose, observation matrix Input matrix ; Given the conjugate of the input matrix ,Depend on Estimate each column The complex values ​​of each column are used to obtain the transmission matrix value at the focal point. .

Images