Stroke monitoring

Abstract

A computer-implemented method for continuously monitoring stroke during a critical recovery period, the method comprising the steps of: (i) accessing initial image data, the initial image data representing an initial image of a subject's brain containing a stroke region; (ii) accessing scattering parameter data, the scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed around the subject's brain; and (iii) processing the scattering parameter data and the initial image data using a gradient-free optimization method to generate an estimate of the spatial dimensions of the stroke region within the subject's brain, wherein the initial image of the subject's brain is used as prior information to improve the accuracy of the generated estimate, and the spatial dimensions of the stroke region are generated as global parameters of the gradient-free optimization method.

Description

Technical Field

[0001] The present invention relates to medical imaging, and in particular to apparatus and methods for continuous monitoring of stroke during the critical recovery period. Background Art

[0002] The World Health Organization defines a stroke as "a neurological deficit of cerebrovascular origin lasting more than 24 hours or interrupted by death within 24 hours". This suggests that rapid action is needed to enable the administration of treatments to reduce the severity of stroke, since the resulting loss of brain function becomes irreversible approximately 6 hours after symptom onset. However, such rapid action is currently being undertaken in partial medical darkness, as there are currently no stroke monitoring tools to effectively measure stroke "geometry" and its response to treatment.

[0003] The first six hours after a stroke occurs are often referred to as the "critical recovery period" ("CRP"), which is defined as the "period during which functional loss caused by the stroke can be recovered." Therefore, once a stroke is diagnosed, medical treatment of the stroke is immediately initiated. However, particularly during the CRP period, medical treatment needs to include essentially continuous monitoring to assess the effects of medications or other treatments on the "geometry" of the stroke and, if necessary, to provide corresponding medical adjustments (e.g., immediate change to a different treatment based on the behavior of the stroke). However, currently available imaging modalities do not currently facilitate such critical needs, for which an example is described in the following literature: P.D. Schellinger et al., "Monitoring intravenous recombinant tissue plasminogen activator thrombolysis for acute ischemic stroke with diffusion and perfusion MRI," Stroke, Vol. 31, No. 6, pp. 1318-1328, March 2000 ("Schellinger").

[0004] Geometric changes of stroke during CRP

[0005] To demonstrate stroke behavior during CRP, Figure 1A A stroke model derived from an MRI 3 hours after symptom onset is shown, similar to the medical situation in Schellinger. The model was derived from images generated by a "stroke MRI" imaging modality used in emergency situations. The main purpose of this imaging modality is simply to reveal the severity of the stroke, so this imaging modality usually has lower specificity than conventional MRI images that require longer computing time to generate.

[0006] Figure 1B and Figure 1C Indicates that 5.5 hours after the onset of symptoms during the CRP Figure 1A Two possible variations of the shape of the stroke. Figure 1B The enlargement of the stroke shown in the Figure 1 is often referred to as "hypoperfusion" and is characterized by insufficient oxygen supply to the brain tissue due to the extension of the "clot", which reduces the dielectric properties of the affected area to values ​​10% to 15% below those of the surrounding tissue. This occurs when thrombolytic therapy (a method of breaking up the "clot" by injecting or directly delivering (using a catheter) clot-dissolving drugs to the affected area) is given too late - usually after 4 and a half hours from symptom onset, and has the opposite effect of causing the stroke area to grow, rather than shrinking as expected.

[0007] on the other hand, Figure 1C Demonstrated the breakdown of clots when thrombolytic therapy is promptly administered, usually within the first 4 hours after stroke onset. Figure 1B and Figure 1C Geometric changes in stroke are key responses to different stroke treatments, but are currently undetectable due to the lack of effective stroke monitoring methods. Figure 1B and Figure 1C geometric changes.

[0008] During a stroke, when a brain artery ruptures, the blood that rushes into the surrounding tissue takes about four minutes to form a clot. Therefore, by the time the patient arrives, the affected area of ​​the brain has already been occupied by a clot.

[0009] Barriers to real-time monitoring

[0010] The lack of continuous monitoring technology stems from some inherent limitations of existing medical imaging modalities, primarily MRI and X-ray. For example, the bulky and static structure of MRI makes it impractical to use as a continuous monitoring tool, as patients with severe stroke are usually in an emergency medical stage, which prevents them from being repeatedly transported between the intensive care unit (ICU) and the imaging unit for monitoring. In addition, the ionizing effect of X-rays and the genotoxic effects of MRI on the human body prevent these imaging modalities from being used to continuously expose the body to potentially unsafe doses of radiation.

[0011] To address these difficulties, electromagnetic tomography (EMT) has been introduced as a possible monitoring mode. EMT involves reconstructing an image of an unknown subject based on measurement data in the form of what is generally referred to as an "S matrix" or "scattering parameters" in the art. This imaging mode benefits from the fact that EMT hardware is lightweight, so it can be carried and installed in every medical unit (including ICU), and this imaging mode also benefits from the fact that EMT hardware is safe because it does not cause the ionizing and genotoxic effects of MRI and X-ray imaging modes.

[0012] However, compared with these mature technologies, EMT has some limitations. In particular, in the UHF (300 MHz to 3 GHz) and S (2 GHz to 4 GHz) frequency bands where medical EMT is usually performed, the sharp corners of human head tissue (especially Figure 1A The presence of highly localized evanescent waves at the rounded corners between the brain and cerebrospinal fluid (shown in Figure 1) and diffraction effects are very significant. Therefore, EMT images are geometrically rough and cannot accurately represent the actual shape of the subject being imaged (e.g., stroke area), which is desirable for CRP monitoring.

[0013] For example, Figure 1D It is shown Figure 1A The relative dielectric constant (ε r ) as generated using eight imaging dipole antennas operating at 1.5 GHz. Figure 1A A rough comparison shows that the computed images do not accurately represent the actual spatial distribution or dielectric properties of the stroke (the retrieved conductivity denoted by σ is excluded from the demonstration because the corresponding retrieval accuracy for conductivity σ is generally lower than that for relative permittivity ε). r The corresponding retrieval accuracy is even worse). In addition, the corresponding computing time required to generate medical images by EMT is quite high. For example, using an Intel Core i7 processor with a CPU of 3.6 GHz TM Computer platform with i7-4790 CPU and 16GB RAM, generated at 15dB signal-to-noise ratio (SNR) Figure 1D The time to acquire an image is about 76 minutes. Therefore, the presence of evanescent waves in the UHF and S bands, diffraction effects, and the prohibitive computational time required preclude the use of conventional EMT for stand-alone stroke monitoring. As proposed by Schellinger: "There remains a need for a rapid stroke imaging tool that is sensitive enough to detect intracerebral hemorrhage (ICH) within the first 6 hours and that, if it exists, can identify tissue at risk."

[0014] It would therefore be desirable to overcome or alleviate one or more difficulties of the prior art or at least provide a useful alternative. Summary of the invention

[0015] According to some embodiments of the present invention, there is provided a computer-implemented method for continuously monitoring a stroke during a critical recovery period, the method comprising the following steps:

[0016] (i) accessing initial image data representing an initial image of a subject's brain including a stroke region;

[0017] (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed about the subject's brain; and

[0018] (iii) processing the scatter parameter data and the initial image data using a gradient-free optimization method to generate an estimate of the spatial dimensions of the stroke region within the subject's brain, wherein the initial image of the subject's brain is used as prior information to improve the accuracy of the generated estimate, and the spatial dimensions of the stroke region are generated as global parameters of the gradient-free optimization method.

[0019] In some embodiments, the spatial dimensions of the stroke region are initially determined by optimizing the spatial dimensions of a first predetermined permittivity value of the stroke region of the subject's brain and the spatial dimensions of a second predetermined permittivity value of a non-stroke region of the subject's brain.

[0020] In some embodiments, the spatial dimensions and relative permittivity of the stroke region are generated as global parameters of a gradient-free optimization method.

[0021] In some embodiments, the shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimensions of the stroke region are determined by determining the spatial dimensions of the overlapping ellipses.

[0022] In some embodiments, the overlapping ellipses have a minor axis with a fixed spatial dimension, and the spatial dimension of the overlapping ellipses is determined as two parameters corresponding to the major axis of the overlapping ellipses. In some other embodiments, the spatial dimension of the stroke region is determined by determining four geometric parameters.

[0023] In some embodiments, the method includes repeating step (ii) and step (iii) over a continuous period of time to monitor the growth or shrinkage of the stroke area over time.

[0024] The non-gradient optimization method may be a Nelder-Mead non-gradient optimization method.

[0025] An initial image of the subject's brain may be generated by magnetic resonance imaging or x-ray imaging or electromagnetic tomography imaging.

[0026] According to some embodiments of the present invention, there is provided a device for continuously monitoring a stroke during a critical recovery period, the device comprising:

[0027] Memory;

[0028] at least one processor; and

[0029] At least one computer-readable storage medium having instructions stored thereon, which, when executed by the at least one processor, cause the at least one processor to perform the following steps:

[0030] (i) accessing initial image data representing an initial image of a subject's brain including a stroke region;

[0031] (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed about the subject's brain; and

[0032] (iii) processing the scattering parameter data and the initial image data to estimate the spatial dimension of the stroke region within the subject's brain, wherein the initial image of the subject's brain is used as prior information to improve the accuracy of the determination, and the spatial dimension of the stroke region is determined as a global parameter of the gradient-free optimization method.

[0033] In some embodiments, the spatial dimensions of the stroke region are initially determined by optimizing the spatial dimensions of a first predetermined permittivity value of the stroke region of the subject's brain and the spatial dimensions of a second predetermined permittivity value of a non-stroke region of the subject's brain.

[0034] In some embodiments, the shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimensions of the stroke region are determined by determining the spatial dimensions of the overlapping ellipses.

[0035] In some implementations, the spatial dimensions of each of the overlapping ellipses are determined as two geometric parameters.

[0036] In some embodiments, the spatial dimensions of the stroke region are determined by determining four geometric parameters.

[0037] In some embodiments, the apparatus comprises repeating step (ii) and step (iii) over a continuous period of time to monitor the growth or shrinkage of the stroke area over time.

[0038] The non-gradient optimization method may be a Nelder-Mead non-gradient optimization method.

[0039] According to some embodiments of the present invention, at least one computer-readable storage medium is provided, on which instructions are stored, and when the instructions are executed by at least one processor of a brain monitoring device, the at least one processor performs the steps of any of the above methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Some embodiments of the invention are described hereinafter, by way of example only, with reference to the accompanying drawings, in which:

[0041] Figure 1A is derived from a cross-sectional plan view MRI image of a subject's head containing the stroke region and shows stroke model information used to initiate EMT, where the corresponding tissue dielectric properties (ε r ,σ) is taken from Gabriel;

[0042] Figure 1B and Figure 1C According to Schellinger Figure 1A The two corresponding possible key responses of stroke during CRP;

[0043] Figure 1D is an image of a subject's head that was generated by a prior art EMT method based on robust gradient-based optimization and was not initiated by a stroke MRI image, and that image clearly cannot accurately identify the stroke region;

[0044] Figure 1E and Figure 1F Is caused by and for Figure 1D The corresponding EMT image was generated in the same way as in Figure 1A The stroke MRI images initiated, Figure 1E Accurately identify the stroke area, Figure 1F Its subsequent response to drug treatment (thrombolysis) is shown at each pixel. However, this improved EMT retrieval method takes more than an hour, making it impractical for stroke monitoring during CRP;

[0045] Figures 2A to 2G The three stroke parameters a0, b0 and ε are shown as r (Geometric parameters a0 and b0 in Figure 1A Different optimization operations of Nelder-Mead (“NM”) optimization of geometric operations in a 3D parameter space (shown in FIG. 4 );

[0046] Figure 3is a flow chart of the NM optimization method according to the described embodiment of the present invention; after constructing the vertex, the method first checks whether the "best" vertex satisfies the truncation condition. If the "best" vertex does not satisfy the truncation condition, a first error reduction operation, namely "reflection", is performed. Therefore, if the reflected vertex has a lower error relative to the best vertex, an expansion operation is performed in the hope of finding an even better vertex. If the reflected vertex does not have a lower error relative to the best vertex, forward / backward tightening and contraction are performed based on the error level of the reflected vertex relative to the worst vertex and the second worst vertex to eventually find a new vertex that is at least better (has a lower error) than the "worst" vertex in each iteration;

[0047] FIG. 4A to FIG. 4C are the corresponding graphs representing the retrieval method and convergence behavior of the NM optimization method described in this paper;

[0048] FIG. 5A to FIG. 5D are corresponding graphs showing the influence of some influencing factors on the accuracy and computation time of NM optimization in stroke monitoring during CRP with 8 antennas and using two stroke shape parameters a0, b0; Figure 5A The effect of the total number of antennas when SNR = 15dB is shown. Figure 5B shows the impact of SNR, Figure 5C The effect of stroke size when SNR = 15 dB is shown. Figure 5D The impact of the difference between each individual tissue and the database in Gabriel is shown when SNR=15dB;

[0049] 6A to 6F : A NM-optimized retrieval method for 3D MRI-derived stroke patterns, where:

[0050] Fig. 6A is a representation of a three-dimensional EMT geometry having a circular antenna array disposed around a patient's head;

[0051] Figure 6B A close-up of a stroke region within a patient's head imaged by the three-dimensional EMT method described herein is shown;

[0052] Figure 6C An imaged stroke area with four geometric parameters defining the three-dimensional stroke area is shown;

[0053] Fig.6D The three-dimensional stroke area imaged during clot breakdown is shown;

[0054] Fig. 6E are plots of the four geometric stroke region parameters and relative permittivity determined by this method as a function of the number of iterations;

[0055] Fig. 6F is a graph showing the convergence of the method (in dB) as a function of the number of iterations;

[0056] Figure 7 is a schematic diagram of a stroke monitoring device according to an embodiment of the present invention; and

[0057] Figure 8 is a flow chart of a stroke monitoring method according to an embodiment of the present invention. DETAILED DESCRIPTION

[0058] The possibility of EMT as a monitoring modality prompted the inventors to develop a new EMT method that is suitable for stroke monitoring during CRP. In fact, although some inherent limitations of EMT mentioned above - namely, the presence of evanescent waves at the UHF band and the S band and diffraction effects - are unavoidable, the inventors determined that the excessively long computational times of the prior art EMT methods are due to the numerical formulation of the gradient-based optimization used in each of the prior art EMT systems. Since these optimizations seek optimal values ​​of "variables" at each pixel of the resulting image, high computational times are generally required to find these optimal values. In other words, if the spatial distribution and dielectric properties of the stroke area in the brain are unknown, they are considered to be variables and their optimal values ​​must be retrieved at each "pixel" of the image (e.g. Figure 1A As shown in ). Since a relatively high number of pixels is required to have a clinically acceptable spatial resolution (even when considering the limited maximum achievable spatial resolution of EMT), the prior art EMT calculation methods are very inefficient for stroke monitoring during CRP. Figure 1D It is obvious that the images computed by these inefficient methods are not accurate enough in any case.

[0059] With this in mind, the inventors realized that since the patient is immediately transferred to the imaging unit upon arrival, it is possible to obtain the initial spatial distribution of the stroke area from an early diagnostic stage using a high-resolution imaging modality such as MRI (or X-ray or even EMT) as prior information before transferring the patient to the ICU (generating an initial image of the stroke area using a stroke MRI imaging modality takes less than an hour), e.g. Figure 1AAlthough MRI (or X-ray or EMT) imaging cannot practically be used for continuous monitoring for the reasons mentioned above (safety, cost, and computational time), stroke MRI (or X-ray or EMT) images generated upon patient arrival can be aligned or "registered" as described in: G. Boverman, CEL Davis, SD Geimer, and PM Meaney, "Image registration for microwave tomography of the breast using priors from nonsimultaneous previous magnetic resonance images," IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, Vol. 2, No. 1, pp. 2-9, March 2018; and R.L. Leijsen, W.M. Brink, C.A.T. van den Berg, AG Webb, and R.F. Remis, "Three-dimensional contrast source inversion-electrical properties tomography," IEEE Transactions on Medical Imaging, vol. 2, no. 1, pp. 2-9, March 2018. Trans. Med. Imag., Vol. 37, No. 9, September 2018 (“Leijsen”), in which Leijsen gave the dielectric properties of different tissues.

[0060] like Figure 1E and Figure 1F As shown, (in this example Figure 1A The registered images can be employed in an otherwise conventional EMT device installed in the ICU (in fact, MRI is performed only once) and then used as a priori information to generate images of stroke response (in Figure 1B and Figure 1C ), showing a significant improvement in retrieval accuracy. The reason for the improvement is that the registered stroke MRI images are used to constrain the space of possible retrievable images generated by the gradient-based optimization, thereby allowing EMT to converge to a more accurate image. However, due to the pixel-based optimization nature of conventional (and therefore gradient-based) EMT, the computational time of this method is still prohibitive, which for this example is approximately one hour. This is significantly lower than the computational time of EMT without MRI images (76 minutes). Figure 1DThis is only a slight reduction in time compared to the example of . Therefore, if gradient-based EMT is initiated by MRI, X-ray or CT scan images, gradient-based EMT may only be used in non-emergency situations, such as tumor screening (as Leijsen mentioned, for the 3D problem, the reported computation time is about 11 hours on a standard computer).

[0061] This article describes a device and method for continuously monitoring cerebral stroke during the critical rehabilitation period (CRP), and for convenience, the device and method are also referred to as stroke monitoring device and method in this article. In the stroke monitoring method and device described in this article, the image of MRI (or X-ray or EMT) registration is implemented as prior information to provide the initial geometry and dielectric properties of the stroke area. However, instead of calculating the inefficient prior art of the variable at each pixel of the imaging area, the shape, dimension and dielectric properties of the stroke area are calculated as global parameters. Therefore, as described below, the global values ​​of these parameters are continuously updated by the gradient-free optimization method at a later monitoring time.

[0062] For example, in the described embodiment, the parameters are two geometric parameters defining the shape and dimensions of the stroke area as the semi-major axes of two overlapping ellipses with a fixed minor axis of 1 cm - e.g. Figure 1A As shown, a0 and b0 (initial values ​​are a0 = 1.5 cm, b0 = 2 cm) - and the relative dielectric constant ε that defines the stroke area r The third parameter of . Considering the stroke shape and dielectric properties as global parameters reduces the number of unknown quantities from e.g. Figure 1A The approximately 1000 pixel-based variables for the image in are reduced to only 3 global parameters. This reformulation enables gradient-free optimization to be used for EMT as an efficient technique for parameter optimization (because the parameters are optimized globally, rather than locally at each pixel). Obviously, this gradient-free optimization cannot provide pixel-based images, such as Figures 1D to 1F However, for medical needs during CRP, rapid and accurate tracking of stroke parameters (e.g., a0, b0, and ε) every two minutes is necessary. r ) is more medically relevant than obtaining a crude picture of stroke on an hourly basis through gradient-based optimization.

[0063] Gradient-free optimization methods were developed to address problems where gradient-based optimization is not applicable, especially when the function to be minimized during optimization is not differentiable or smooth.

[0064] Various gradient-free optimization methods have been developed for different electromagnetic and antenna applications, including Nelder-Mead ("NM"), genetic algorithms ("GA"), and particle swarm optimization ("PSO") methods. In cases where computational time is the primary consideration in the optimization method in addition to accuracy, the NM optimization method is typically the fastest gradient-free optimization method. In the case of stroke monitoring during CRP, time is life, and the inventors believe that NM optimization is most suitable for monitoring requirements (although other gradient-free optimization methods may be used in other embodiments). Since NM optimization has not been used for EMT applications before, an overview of the NM optimization method for this particular application is described below. General NM methods are described in the following literature: N. Pham, A. Malinowski, and T. Bartczak, "Comparative study of derivative free optimization algorithms", IEEE Trans. Industr. Inform., Vol. 7, No. 4, pp. 592-600, November 2011 ("Pham").

[0065] In the described embodiment, Figure 7 As shown, the described method is performed by a stroke monitoring device. In use, the device communicates with a vector network analyzer (VNA) or transceiver 701, which in turn is connected to a microwave antenna array 705.

[0066] As shown, the microwave antenna array 705 is arranged to receive the head 704 of a patient whose brain is to be imaged, so that each antenna in the array can be selectively excited to radiate electromagnetic waves or microwave frequency signals into the subject's head and through the subject's head to be scattered, and the corresponding scattered signals are detected by all antenna arrays (including the antenna that transmits the corresponding signal).

[0067] It will be apparent to those skilled in the art that a vector network analyzer (VNA) 701 excites the antennas as described above and records the corresponding signals from the antennas as data (referred to herein as "scattering" data) representing the amplitude and phase of the scattered microwaves in the form of what are referred to in the art as "scattering parameters" or "S parameters". The VNA 701 sends this data to a device for processing to generate information about the internal features of the imaged object (e.g., brain clots, bleeding sites, and other features). In the described embodiment, a VNA with a large dynamic range of greater than 700 dB and a noise floor of less than -700 dBm can be used to activate the antennas to transmit electromagnetic signals over a frequency band of 0.5 GHz to 4 GHz and receive scattered signals from those antennas.

[0068] Although the apparatus of the described embodiments is in the form of a computer, this need not be the case in other embodiments. Figure 7 As shown, the stroke monitoring device of the described embodiment is a 64-bit Intel architecture computer system, and the stroke monitoring method performed by the stroke monitoring device is implemented as programming instructions of one or more software modules 702 stored on a non-volatile (e.g., hard disk or solid state drive) storage device 704 associated with the computer system. However, it is obvious that at least a portion of these methods may alternatively be implemented in one or more other forms, such as being implemented as configuration data of a field programmable gate array (FPGA), or being implemented as one or more dedicated hardware components such as an application specific integrated circuit (ASIC), or being implemented as any combination of such forms.

[0069] The stroke monitoring device includes a random access memory (RAM) 706, at least one processor 708, and external interfaces 710, 712, 713, 714, all of which are interconnected by a bus 716. The external interfaces include: a network interface connector (NIC) 712, which connects the stroke monitoring device to a communication network such as the Internet 720; a universal serial bus (USB) interface 710, at least one of which can be connected to a keyboard 718 and a pointing device such as a mouse 719; and a display adapter 714, which can be connected to a display device such as an LCD panel display 722.

[0070] The stroke monitoring device also includes an operating system 724 such as Linux or Microsoft Windows, and in some embodiments includes additional software modules 726 to 730, including: web server software 726, such as Apache, which can be run on http: / / www.apache.org scripting language support 728, such as PHP, available at http: / / www.php.net, or Microsoft ASP; and structured query language (SQL) support 730, such as MySQL, available from http: / / www.mysql.com, which allows data to be stored in and retrieved from a SQL database 732.

[0071] The web server 726, scripting language module 728, and SQL module 730 together provide the stroke monitoring device with general capabilities to allow remote users with standard computing devices equipped with standard web browser software to access the stroke monitoring device, particularly to monitor the progression of stroke during CRP. A. NM Optimization of EMT: Initialization

[0072] for Figure 1B and Figure 1C (Depend on Figure 1A In order to solve the EMT problem shown in FIG. 1 , eight imaging dipole antennas are arranged in a circular configuration around the subject's head so that a 2D MRI model of the subject's head passes through the ports of the dipole antennas (because the S matrix implemented in the optimization method is recorded at the dipole ports).

[0073] The antenna irradiates the head with a sinusoidal electromagnetic wave at 1.5 GHz, say 5 and a half hours after symptom onset, at which time thrombolysis has been administered. This timing example is taken from the stroke case described by Schellinger, where the subject's arrival time was approximately 3 hours after symptom onset and a 2D single-slice stroke MRI was prepared 0.75 hours after arrival. The antenna then records the corresponding scattered fields in the form of an S matrix (the effect of the number of antennas on accuracy and retrieval time is discussed below). This S matrix is ​​then implemented to update the three stroke parameters a0, b0, ε r , the three stroke parameters are Figure 1A The initial value in is X0 = (a0 = 1.5 cm, b0 = 2 cm, ε r =39), which can be considered to define points (or vertices as described below) in 3D space (since the number of parameters is 3).

[0074] These parameters were chosen for the following reasons. Since the relative permittivity ε of the stroke region through the subject's brain r changes to reflect the effect of any drug administered to the stroke patient, so the relative dielectric constant ε r can be considered as a reliable parameter for monitoring stroke and converges to a value matching the S matrix at a later time. As mentioned above, due to the relative permittivity ε r In comparison, the retrieval accuracy of the conductivity σ of the stroke area is usually poor, so the conductivity is excluded from the retrieval. In addition, the geometric parameters or shape parameters a0, b0 are selected to best reflect the geometric changes (extension or contraction) of the stroke area. Figure 1A The smooth stroke shape shown, a0, b0 are the semi-major axes of two ellipses, which are considered to most accurately represent the geometric changes in the stroke area relative to possible alternative geometric parameters (more complex geometries are discussed below). The inventors have found that these shape parameters also continuously converge to the values ​​that best fit the corresponding S matrix. Figure 1A The 2D physical problem domain in Figure 2A The 3D mathematical space of the parameters constructed in is different.

[0075] In addition to X0, the range of variation of the three parameters defines the available parameter space. Assuming that the stroke area suffering from hypoperfusion ( Figure 1B) can extend to the entire left hemisphere of the brain, and successful clot resolution ( Figure 1C ) can completely remove the clot (i.e., shrink the stroke area to zero). This large range corresponds to 0≤a0≤4cm, 0≤b0≤7cm, and 39≤ε r ≤43.5 in the parameter space, indicating that even without reliable predictions about stroke response to drugs, the optimization method does not require a very precise parameter range when searching for the optimal vertex. Figure 2A As shown, in order to optimize these three parameters, NM optimization first constructs an isometric simplex (generalized triangle) with 3+1 vertices in the 3D space of parameters. In fact, one of the vertices is the initial guess (i.e., starting point) X0 provided by the initial MRI image. All other vertices with equal distances c are derived as follows by adding the following vectors to this initial guess:

[0076]

[0077] The vector components are:

[0078]

[0079] and where N is the number of parameters (i.e., 3 in the described embodiment). Typically, c=1 allows the method to search in a sufficiently large volume at the initial step. Small values ​​of c typically require long computation times and may mislead the method into finding only local minima. Furthermore, since the presence of noise may cause slightly different vertices to produce the same frequency response (S matrix), having the vertices located far enough apart from each other (c ≥ 1 apart) makes the method robust to noise at early iterations.

[0080] With the simplex constructed, the next step is to evaluate the objective function (i.e., the function to be minimized by NM optimization) at each vertex X0, X1, X2, X3. Figure 1B or Figure 1C In the EMT problem of , the objective function is defined as the L2 norm mismatch between the measured (“meas.”) S matrix and the retrieved (“retr.”) S matrix as follows:

[0081]

[0082] Where Ω represents the imaging domain (in cylindrical coordinates ), χ is defined as Comparison of the dielectric properties of a human head, where is the complex dielectric constant, is the free space dielectric constant, and ω = 2π × 1.5 GHz is the angular frequency. In the described example with eight antennas, the size of the S matrix is ​​therefore 8 × 8; the mismatch is therefore S meas .With S retr The difference between the corresponding matrix elements of . As described in the following literature, E is Figure 1B or Figure 1C The total electric field over the imaging domain of: Chapter 9 of W.C. Chew, “Waves and Fields in Inhomogeneous Media,” IEEE Press, New York, 1995 (“Chew”) and A. Afsari, A. Abbosh, and Y. Rahmat-Samii, “Modified Born iterative method in medical electromagnetic tomography using magnetic field fluctuation contrast source operator,” IEEE Transactions on Microwave Theory and Techniques, DOI: 10.1109 / TMTT.2018.2876228, and E is given by:

[0083]

[0084] Among them, E inc is the incident electric field when there is no object in Ω. Finally, G is the dyadic Green's function given in Chapter 1 of Chew, and is the free space wave number, where μ0 = 4π × 10 -7 is the free space permeability.

[0085] like Figure 2B As shown, after evaluating the objective function of all vertices, three vertices have special importance in the NM technique. The vertex that gives the minimum value of the objective function (given by equation (3) in the described embodiment) is called the "best" vertex and is denoted by X b The vertex with the highest L2 norm mismatch is called the “worst” vertex X w In order to steer the method in a direction that always reduces the mismatch in equation (3), the method determines the auxiliary point with the second highest value of equation (3). This vertex is called the "second worst" vertex X sw Then divide by X w The average value of all vertices except a .X a With X w The line segment between (LXwXa ) is always downward and contains some useful points where equation (3) may have a higher value than X during the optimization method. w Lower mismatch. According to this downward direction search for the optimal point, this method is also called "downhill" optimization.

[0086] After the above initialization, the method performs at least two and at most five "error reduction" operations. In each of these steps, X w The old value of is removed from computer memory (i.e., not stored for the next iteration), and all other vertices are rearranged to provide X b , X w , X sw The following operations are described below:

[0087] Reflection: The first optimization step in the NM method is to reflect the worst vertex X on L with the same length as follows:

[0088] X r =2X a -X w (5)

[0089] like Figure 2C This operation is to check Whether the movement in the direction should be continued by the method or the other direction to guide the method to the best vertex. If equation (3) is in X r Relative to X b Has a lower mismatch, that is, F(X r )<F(X b ), then before replacing:

[0090]

[0091] This method evaluates the chance of finding an even better vertex (where the programming convention A←B means that the old value A is replaced by the new value B). To do this, always Move further up to perform the expansion operation.

[0092] Extension: According to Figure 2D , the method follows Move further with the same step length, i.e.

[0093] X e =2X r -X a (7)

[0094] Then, equation (3) is also evaluated at this extension vertex. If its value is lower than X b (Even if it is better than X r Worse), that is, F(Xe )<F(X b ), the method performs the following substitution:

[0095]

[0096] and (iteratively) return to the reflection step. The method does not immediately accept X even though it is the best vertex found among the other vertices. b ←X r The reason for this comes from the fact that since the vertex is located at the e The interior of the new simplex formed by X is therefore retained by this method. Therefore, by performing the expansion, X r The neighborhood of is only protected as a subdomain, where there may be some other good vertices or even better vertices to minimize equation (3). However, if equation (3) is e No less than X b , then the substitution of equation (6) is performed and the method iteratively returns to the first operation (ie, reflection).

[0097] Forward contraction: Equation (6) or Equation (8) assumes that F(X r )<F(X b ). If this is not implemented, but F(X sw )<F(Xr)<F(X w ), then this method has been direction is shifted too much, and a better vertex may be located at a position other than X r Therefore, if Figure 2E As shown, by r To X a Returns half a step size (usually a balanced step size) to execute along The forward contraction, that is

[0098] X fc =1.5X a -0.5X w (9)

[0099] If F(X fc )<F(X r ), then by returning to the initialization step and rearranging the vertices from the worst vertex to the best vertex, at vertices X, X b , X sw , X fc A new simplex is formed on it.

[0100] Tighten backwards: Figure 2F As shown, if F(X w )<F(X r ),but There may still be some vertices that can improve the mismatch of equation (3), and as follows, in the same way as forward contraction, by a To X w Perform a backwards contraction by moving half a step backward:

[0101] X bc =0.5X a +0.5X w (10)

[0102] If F(X bc )<F(Xw ) , then by returning to the initialization step and rearranging the vertices from worst to best, b , X sw , X bc The above structure is Figure 2F The new simplex is shown.

[0103] Contraction: Nevertheless, if none of the above conditions occur, the final step to find a better direction towards the best vertex is to contract the simplex. To do this, keep only the best vertex X b , and for the other vertices, the contraction operation is performed as follows (for each i-th vertex):

[0104] X i (new) = 0.5X b +0.5X i (old) (11)

[0105] The method then returns to the initialization step to rearrange the Figure 2G The new vertices formed in the contraction step are shown. This iterative process of NM optimization is continued until F(X b ) satisfies the truncation condition of the objective function, which in the described embodiment is:

[0106] F(X b )<10 -7 (12)

[0107] X that satisfies the cutoff condition b The value of is stored as the final result. The cutoff condition in equation (12) is chosen to be very small to ensure the accuracy required to meet the retrieval parameters. Larger values ​​of the cutoff condition do not result in very accurate parameter retrieval. To demonstrate all these steps at once, Figure 3 is a flow chart of the NM optimization method for EMT, and the corresponding pseudo code is given below.

[0108]

[0109]

[0110] I. Practical NM Gradient-free Optimization: 2D Retrieval

[0111] To effectively monitor what is known as hypoperfusion and Figure 1B Different stroke behaviors shown in Figure 1C The above NM optimization method is used to optimize the three parameters a0, b0 and ε r . Only shape parameters a0 and b0 are selected for optimization and excluding Figure 1B or Figure 1C The reason for the position parameter of the coordinates of the center of the stroke comes from the fact that during hypoperfusion or clot breakdown, the stroke position does not change. Instead, the stroke only expands or shrinks around its center.

[0112] A.NM Optimization Performance

[0113] For low perfusion, practical values ​​of the parameters in the given example are a0 = 2.5 cm, b0 = 3 cm, ε r =39, and the retrieved value - in Figure 4A Its evolution in successive iterations is shown in r =39.08. Figure 5A As shown, the NM method gradually becomes stable as the simplex becomes smaller during the iteration. Figure 1C The actual and retrieved values ​​of the NM optimization method are a0 = 1 cm, b0 = 1 cm, ε r =40 and a0 = 0.997 cm, b0 = 1.002 cm, ε r =40.07( Figure 4B ). It should be noted that during clot dissolution, the relative permittivity of the affected area approaches that of the brain again upon recovery. In both cases, despite the high noise of the simulated environment (the SNR of the S matrix was set to 15 dB), the parameters were still accurately updated. To provide a quantitative analysis of this accuracy, the “retrieval error” was defined as:

[0114]

[0115] In each case, the retrieval error was low: -18.86 dB for the low perfusion result and -27.95 dB for the clot breakdown result.

[0116] However, the key point is the retrieval time. Figure 4C, once the termination condition of equation (12) is met, the method terminates. In the described example, the hypoperfusion and clot decomposition retrievals required 94 and 69 iterations, respectively, and each iteration took 1.25 seconds. Thus, within every two-minute time frame, the stroke is continuously monitored to evaluate treatment methods and to prepare for craniotomy (surgery to remove a portion of the subject's skull to directly access the affected area) immediately if needed.

[0117] Figure 4A and Figure 4B A graphical illustration of the performance of the NM method is provided for only the first 70 iterations and the first 40 iterations, respectively (the changes in the parameters retrieved in further iterations are too small to be usefully plotted in the same figure). Figure 1A ) results in higher retrieval errors and longer convergence times compared to clot decomposition retrieval.

[0118] A. Factors affecting accuracy and computation time

[0119] Among the different factors that affect the accuracy (retrieval error) and computation time of NM-optimized EMT, the main contributors are: the total number of imaging antennas, SNR, shape parameters a0, b0, and slight differences in the dielectric properties of each subject's head tissue relative to the database described in the following literature, which are used to align the initial stroke MRI images; the literature is as follows: C. Gabriel, S. Gabriel and E. Corthout, "The dielectric properties of biological tissues: I. Literature survey", Phys. Med. Biol., Vol. 41, No. 1, pp. 2231-2249, 1996 ("Gabriel"). To describe these effects, since the retrieval of low perfusion is more challenging (with higher retrieval error), the impact of the number of antennas, SNR and Gabriel's accuracy level is described for this medical scenario.

[0120] Figure 5AIt is shown that increasing the number of antennas to a certain level significantly improves the accuracy at the expense of increasing the computation time for each iteration. This improvement comes from the fact that having additional antennas is equivalent to having more information. The total number of iterations is reduced from 94 to 78 for the same reason, because the space of possible solutions is now more restricted and the method requires fewer iterations to satisfy the truncation condition. However, the computation time for each iteration increases dramatically due to the increase in the size of the objective function, which corresponds to the size of the S matrix. In practice, having M imaging antennas results in an M×MS matrix. In the literature, 8 are suggested as the minimum number of antennas to provide medical images of the head with acceptable accuracy.

[0121] The second factor to consider is SNR. When the measured data is highly contaminated by noise, vertices that are close to each other (in Figure 2A have close parameter values ​​in the same frequency response that is most likely dominated by noise. Therefore, the accuracy of this method will decrease because it can no longer distinguish between close vertices. Figure 5B As shown, if the domain of medical imaging is highly isolated from noise, both the retrieval accuracy and the number of required iterations are improved. This is because the measured S matrix in equation (3) contains more reliable information for retrieving stroke responses (hypoperfusion), so the vertex that best minimizes the objective function in equation (3) is closer to the actual vertex. The computation time for each iteration remains unchanged, i.e., 1.25 seconds, because it is a function of the number of imaging antennas and the retrieval parameters.

[0122] The next factor to consider is the size range of the stroke region that can be accurately retrieved. To this end, based on experience, a cross-sectional factor a0×b0 is introduced to represent the stroke size. Figure 5C The following is shown: If a0×b0<0.64cm 2 Or a0×b0>9cm 2 , then for Figure 1B As mentioned above, the diffraction effects and the presence of evanescent waves are significant in the UHF band and S band, and these effects become extremely destructive when the size of the stroke area becomes smaller than a quarter wavelength. When the wavelength at 1.5 GHz is 3.2 cm inside the clot (where ε r =39), a quarter wavelength is 0.8 cm, and we get 0.64 cm 2 On the other hand, for very high section factors, for example, a0×b0>9cm 2, the change of stroke shape cannot be accurately tracked because it inevitably requires additional shape parameters (besides a0, b0) to more accurately model the irregularity and asymmetry of the stroke, especially when the stroke deforms from one side through the skull but can still extend from the other side. From the perspective of computational time, if only two shape parameters are used for large-sized stroke regions, the method is not suitable for Figure 1B The stroke in requires more iterations to make the initial shape of the stroke partially conform to the extended shape, e.g. Figure 5C This is due to the fact that adding each shape parameter and its corresponding range of variation will further restrict the space of possible solutions. Treating only two shape parameters does require more iterations to partially match the irregularity and asymmetry of the extended stroke (this is explained below). However, it is important to note that increasing the number of stroke parameters as a solution to improving retrieval error and reducing total iterations does not necessarily reduce computational time, because adding each parameter to the problem adds another dimension to the parameter space. Therefore, the computational time for each iteration increases exponentially.

[0123] The last factor is the influence of the level of consistency between the database in Gabriel and the dielectric properties of each individual person's head. Since the gradient-free optimization method focuses only on retrieving global stroke parameters, it is very important to ensure that the dielectric properties in Gabriel (through which the stroke MRI is registered to the EMT device) are highly accurate to avoid significant retrieval deviations from the expected values. Fortunately, this consistency is usually very high due to the construction Figure 1A The material of the tissues in is identical for different individuals and the small differences between these materials and the database in Gabriel are mainly due to partial statistical randomness. Taking this into account, Figure 5D right Figure 1A A random difference of up to 5% was applied to each tissue in - relative to the corresponding value for each tissue reported in Gabriel - to demonstrate the robustness of the method to small differences in the dielectric properties of different head tissues. As can be seen, both the computational speed and the retrieval error decreased slightly, but still fell within the acceptable requirements of CRP.

[0124] By studying the influence of influencing factors on the accuracy and computational time of the NM optimization method in medical EMT applications, the method is applied to a more complex problem in the following, where, according to Leijsen FIG. 4A to FIG. 4C , 3D MRI-derived stroke models complement NM methods, such as 6A to 6F Four shape parameters are implemented in 3D to monitor more complex clot disintegration processes. I. NM gradient-free optimization in practice: 3D retrieval

[0125] For real 3D EMT problems, the accuracy of the retrieved parameters can be further improved if well-developed numerical methods such as finite element modeling (FEM) are used to simulate the total electric field in the domain or the retrieved S-matrix. This can be achieved by directly solving the wave equation and considering the entire 3D physical structure of the imaging antenna (see Fig. 6A ) rather than using the point source based Equations (3) and (4) given above. In fact, the reason for implementing the point source based objective function of Equation (3) given above is that in EMT, the same objective function is used in the prior art gradient based optimization methods. Therefore, such a point source approximation provides a better basis for comparing the two methods. However, due to 6A to 6F The problem shown has the additional complexity of a real antenna structure and is simulated using the commercial COMSOL software package. Figure 6B As shown in Figure 2, the stroke region within the entire head model was identified as a complex shape defined by triangular mesh elements, and the stroke parameters were Figure 6C is shown as having an initial value X′0=(a1=0.41 cm, b1=0.4 cm, c1=0.41 cm, d1=0.42 cm, ε r =39), which are the distances from the center of the stroke area to the four farthest points on the perimeter of the area.

[0126] As in Fig.6D As seen in Figure 2, when thrombolytic therapy is given in time, the clot disintegration method begins. This method will slightly or significantly change the stroke parameters to X1 = (a1 = 0.4 cm, b1 = 0.395 cm, c1 = 0.39 cm, d1 = 0.41 cm, ε r =42).

[0127] according to Fig. 6E The shape parameters of the retrieved stroke area are X1 = (a1 = 0.401 cm, b1 = 0.395 cm, c1 = 0.392 cm, d1 = 0.411 cm) and the dielectric constant of the retrieved stroke area is ε r =41.43, which is very consistent with their actual values. Fig.6D As shown, the changes in the spatial dimensions of the stroke area due to clot disintegration are largest for shape parameters c1, d1, which is why they require additional iterations to approach their actual values. Fig. 6Fis the convergence plot of the NM optimization method. With sufficient shape parameters, less than 50 iterations were required to reach the cutoff condition. In COMSOL Multiphysics installed on the same PC, each iteration took 15 seconds to complete due to the large size of the 3D problem domain. In summary, the behavior of the stroke can be monitored every 11 minutes (i.e., the described device can generate a Fig. 6E A new set of graphs (a), which is a reasonable time for 3D monitoring during CRP. Taking the results described in Pham as ground truth, when using GA or PSO gradient-free methods instead of NM methods, at least 10 times more computational time is required, for example, for this particular problem, the computational time is about 110 minutes. This time is close to the computational time of prior art gradient-based optimization methods and generally exceeds the CRP time range. Therefore, the inventors believe that, in gradient-free methods and general-purpose computer hardware available at the time of writing, an EMT system based solely on NM gradient-free optimization can be effectively used as a 2D or 3D stroke monitoring tool during CRP.

[0128] The EMT monitoring method and device described herein and based on Nelder-Mead gradient-free optimization provide the ability to monitor the expansion or contraction of stroke during CRP, and can therefore potentially be considered as a translational medical advance that increases the chance of survival of stroke. The results described herein show that the method is very efficient in retrieving 2D stroke responses within every 2 minutes or 3D stroke responses within every 11 minutes on a general-purpose computer platform, while other gradient-free methods such as GA or PSO can produce the same output, but at the cost of longer computing time. The described method can be initiated by stroke MRI data obtained from early diagnosis. Then, the shape and dielectric properties (real part of dielectric constant) of the stroke area are defined as global parameters. Subsequently, a small number of imaging antennas can be used to continuously image the patient through a portable EMT system as described herein, and the global parameters are updated using the S matrix recorded by these antennas at each imaging step, thereby identifying the expansion or contraction of the stroke area, particularly in response to one or more treatments. Therefore, the described method and device can improve the treatment process, and therefore improve the survival chance of stroke victims.

[0129] Many modifications will be apparent to those skilled in the art without departing from the scope of the invention.

Claims

1. A computer-implemented method for continuous monitoring of stroke during a critical recovery period, the method comprising the following steps: (i) accessing initial image data representing an initial image of a subject's brain including a stroke region; (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed about the subject's brain; as well as (iii) processing the scattering parameter data and the initial image data using a gradient-free optimization method to generate an estimate of the spatial dimensions of a stroke region within the subject's brain, wherein the initial image of the subject's brain is used as prior information to improve the accuracy of the generated estimate and the spatial dimensions of the stroke region are generated as global parameters of the gradient-free optimization method.

2. The method according to claim 1, wherein: The spatial dimensions of the stroke region of the subject's brain are initially determined by optimizing the spatial dimensions of a first predetermined permittivity value of the stroke region and the spatial dimensions of a second predetermined permittivity value of the non-stroke region of the subject's brain.

3. The method according to claim 1 or 2, wherein: The shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimension of the stroke region is determined by determining the spatial dimension of the overlapping ellipses.

4. The method according to claim 3, wherein: The overlapping ellipses have a minor axis having a fixed spatial dimension, and the spatial dimension of the overlapping ellipses is determined as two parameters corresponding to the major axis of the overlapping ellipses.

5. The method according to claim 1 or 2, wherein: The spatial dimensions of the stroke region are determined by determining four geometric parameters.

6. The method according to claim 1 or 2, comprising repeating step (ii) and step (iii) over a continuous period of time to monitor the growth or shrinkage of the stroke area over time.

7. The method according to claim 1 or 2, wherein: The non-gradient optimization method is the Nelder-Mead non-gradient optimization method.

8. The method according to claim 1 or 2, wherein: The spatial dimensions and relative permittivity of the stroke region are generated as global parameters of the gradient-free optimization method.

9. The method according to claim 1 or 2, wherein: An initial image of the subject's brain is generated by magnetic resonance imaging or x-ray imaging or electromagnetic tomography imaging.

10. A device for continuous monitoring of stroke during a critical recovery period, the device comprising: Memory; at least one processor; as well as At least one computer-readable storage medium having instructions stored thereon, which, when executed by the at least one processor, cause the at least one processor to perform the following steps: (i) accessing initial image data representing an initial image of a subject's brain including a stroke region; (ii) accessing scattering parameter data representing microwaves scattered by the subject's brain and originating from a plurality of antennas disposed about the subject's brain; as well as (iii) processing the scattering parameter data and the initial image data to estimate the spatial dimensions of a stroke region within the subject's brain, wherein the initial image of the subject's brain is used as prior information to improve the accuracy of the determination, and the spatial dimensions of the stroke region are determined as global parameters of a gradient-free optimization method.

11. The device according to claim 10, wherein: The spatial dimensions of the stroke region of the subject's brain are initially determined by optimizing the spatial dimensions of a first predetermined permittivity value of the stroke region and the spatial dimensions of a second predetermined permittivity value of the non-stroke region of the subject's brain.

12. The device according to claim 10 or 11, wherein: The shape of the stroke region is approximated by overlapping ellipses in a two-dimensional plane, and the spatial dimension of the stroke region is determined by determining the spatial dimension of the overlapping ellipses.

13. The device according to claim 12, wherein: The spatial dimensions of each of the overlapping ellipses are determined as two geometric parameters.

14. The device according to claim 10 or 11, wherein: The spatial dimensions of the stroke region are determined by determining four geometric parameters.

15. The device according to claim 10 or 11, wherein: The instructions cause the at least one processor to repeat steps (ii) and (iii) over a continuous period of time to monitor growth or shrinkage of the stroke region over time.

16. The device according to claim 10 or 11, wherein: The non-gradient optimization method is the Nelder-Mead non-gradient optimization method.

17. The device according to claim 10 or 11, wherein: The spatial dimensions and relative permittivity of the stroke region are generated as global parameters of the gradient-free optimization method.

18. The device according to claim 10 or 11, wherein: An initial image of the subject's brain is generated by magnetic resonance imaging or x-ray imaging or electromagnetic tomography imaging.

19. At least one computer-readable storage medium having instructions stored thereon, which, when executed by at least one processor of an apparatus for continuously monitoring stroke during a critical recovery period according to any one of claims 10 to 18, causes the at least one processor to perform the steps of a method according to any one of claims 1 to 9.

Images