Abstract
This paper deals with the possibility of adopting microwave imaging to continuously monitor a patient after the onset of a brain stroke, with the aim to follow the evolution of the disease, promptly counteract its uncontrolled growth, and possibly support decisions in the clinical treatment. In such a framework, the assessed techniques for brain stroke diagnosis are indeed not suitable to pursue this goal. Conversely, microwave imaging can provide a diagnostic tool able to follow up the disease’s evolution, while relying on a relatively low cost and portable apparatus. The proposed imaging procedure is based on a differential approach which requires the processing of scattered field data measured at different time instants. By means of a numerical analysis dealing with synthetic data generated for realistic anthropomorphic phantoms, we address some crucial issues for the method’s effectiveness. In particular, we discuss the role of patient-specific information and the effect of inaccuracies in the measurement procedure, such as an incorrect positioning of the probes between two different examinations. The observed results show that the proposed technique is indeed feasible, even when a simple, nonspecific model of the head is exploited and is robust against the above mentioned inaccuracies.
1. Introduction
A brain stroke occurs when cerebral blood circulation fails as a consequence of a blocked or burst blood vessel, causing a ischemia or an hemorrhage, respectively. Currently, stroke is the second cause of mortality worldwide and a major cause of permanent disability [1]. Moreover, because of the ageing population, such a burden is expected to increase in the future, with a significant, also economical, impact on health-care systems [2]. In such a framework, the most assessed diagnostic tools are magnetic resonance imaging (MRI) and computerized tomography (CT), which provide highly reliable images and indeed play a primary role in stroke management protocols [3].
In the recent years, with the aim to further enhance recovery rate and reduce consequences of strokes, continuous postevent monitoring of physiological parameters in the acute stage has gained an increasing importance [4, 5]. As a matter of fact, recovery and treatment depend on close clinical observation, especially during the first few hours after the onset of stroke [6]. However, CT and MRI cannot contribute to pursue this goal, because they are time consuming, not portable, and cost-ineffective. Moreover, CT uses ionizing radiations, which are harmful for the patient's health and hence not suitable for continuous monitoring. For these reasons, there is an interest in developing cooperating techniques, capable of overcoming these limitations and increase the effectiveness of medical treatments.
Thanks to the evidence that an electric contrast exists between healthy and stroke tissues [7], microwave imaging (MWI) can be exploited in this framework. In addition, microwaves are nonionizing and the involved instrumentation is relatively low cost and portable if compared to CT or MRI. Hence, MWI devices can be safely adopted for repeated examinations and even deployed at the bed of the patient. Therefore, MWI may represent a valid complementary diagnostic tool for the follow-up stage of brain stroke. However, retrieving the map of the electromagnetic properties of brain tissues from the data gather by a MWI device requires to tackle a nonlinear and ill-posed inverse scattering problem [8]. As well known, solution procedures for this class of problems are computationally demanding, thus being not suitable to provide the needed real-time or quasi-real-time results. Moreover, they are affected by the occurrence of false solutions [9], that is, imaging results completely different from the ground truth that impair the reliability of the diagnosis.
The specific nature of the imaging task to be faced in stroke follow-up suggests that the above problems can be overcome. As a matter of fact, the actual target of such a monitoring is the (typically small) variation of the stroke region extent, with respect to the one imaged at the first diagnosis. From the point of view of the MWI processing, this circumstance is indeed a beneficial one, as it entails that the inverse problem to be solved can be well approximated with a linear one, which is free from false solutions and efficiently handled using reliable and assessed inversion strategies.
In this paper, we consider a differential imaging approach for brain stroke followup. In particular, the method provides an estimate of possible variations of the stroke region by processing the difference between monochromatic scattered field data collected at different time instants after stroke diagnosis. In particular, we adopt a Truncated Singular Value Decomposition (TSVD) scheme [10] and, motivated by the encouraging preliminary results [11], we investigate its performance against a priori information on the patient under examination and the presence of inaccuracies in the measurement process. As far as the first point is concerned, the aim is to assess the capability of the proposed method of monitoring the disease's evolution with minimal requirements in terms of patient-specific information. As a matter of fact, while such information could be obtained by converting a high resolution MRI or CT image of the patient's head into relative permittivity and conductivity maps, such a translation may be a difficult task, because of the need of coregistering the images and due to the variability of electric properties of human tissues from patient to patient. Moreover, minimizing patient-specific information allows image processing to be carried out in real time. As far as the second point is concerned, we have analyzed the effect of the presence of uncertainties on the measured data as well as of mismatches occurring in the probes positioning between two measurement stages. As a matter of fact, in follow-up applications data are gathered at different time instants and processed assuming that the overall conditions hold unchanged, so that it is important to consider the effect of such inaccuracies.
It is worth recalling that the concept of differential MWI for medical applications is not new in the literature [12, 13], and some example exists also in the case of brain stroke imaging [14, 15]. However, the approach presented herein is the first one that is explicitly meant for stroke followup. In addition, previous papers concerned with MWI for brain monitoring [8, 14, 15] assume the availability of an almost exact, patient-specific knowledge of the actual scenario before the stroke onset, which poses the above mentioned issues. Finally, as far as application of microwaves for brain stroke classification is concerned, the efforts of Fhager and colleagues which aimed at developing a system using microwave signals to discriminate between ischemic or hemorrhagic events are worth to be mentioned [16–18]. As a matter of fact, the kind of system they have developed, even though meant for a purpose different than imaging, is consistent with the concept we are presenting in this paper, so that it represents in some sense a starting point for the development of a MWI device for brain stroke followup.
The paper is structured as follows. In Section 2, the mathematical background is given and the linear inverse scattering problem at hand is stated. In Section 3, the adopted TSVD inversion procedure is described, with a detail on the role of patient specific information. Section 4 reports a numerical analysis carried out against synthetic data generated for a 2D anthropomorphic head model, concerning the performance of the proposed method against uncertainties on the reference model and on the gathered data. In Section 5, a 3D validation with synthetic data of the overall concept is presented. Conclusions are drawn in Section 6. Throughout the paper a time factor is assumed and omitted.
2. Mathematical Background and Statement of the Problem
Let us suppose that a patient has to be monitored over time after the stroke onset by means of a MWI device. To this end, the scattered field data (encoding the variation with respect to an assumed reference scenario) is measured at different time instants. In particular, the data are collected by means of an array of antennas located on a surface surrounding the head and we assume that a single frequency multiview-multistatic measurement configuration is exploited. In such an arrangement, each antenna acts as both transmitter and receiver, and the field induced when the th antenna is active is recorded at all the probes positions. Notably, owing to reciprocity, it is not necessary to carry out measurements, since only are actually independent and therefore needed to build a complete data-set.
The choice of considering single frequency data is motivated by two reasons. First, our previous studies revealed that the useful bandwidth for microwave brain imaging applications is actually narrow [14], suggesting that processing of multifrequency data would not considerably improve the achievable results, as indeed indirectly confirmed by the fact that numerical simulations in [8] use monochromatic data and by the narrow frequency operation of the experimental device designed by Fhager and coworkers [18]. The second not less important motivation is that monochromatic data allows us to avoid the problem of modeling the frequency dependent electric behavior of tissues and to keep the complexity of the experimental apparatus as low as possible, the latter being an important point to speed-up data acquisition and henceforth the overall examination time.
Under the above assumptions, the time-harmonic scattered field measured at the generic time instant by the receivers located in when the probe in is radiating can be expressed as where represents the investigated region, represents the Green's function for the assumed reference scenario, and is the total field in the reference scenario at the considered time instant. The function is the electric contrast and accounts for the variation at the considered time instant of the electric properties of the imaged scenario, with respect to those of the reference one. For instance, assuming that the reference scenario is the healthy brain and that the instant follows the stroke onset, the contrast represents the stroked area that is the portion of the brain wherein a change in the electric properties has occurred due to the presence of injured tissues.
To express the variation of the stroke region's extent, which is the actual target of followup, let us suppose to have gathered data at two different time instants: at the time of the first diagnosis, say , and after a certain time, say . Owing to the linearity of Maxwell's equations, the resulting differential field is given by
Such an expression is a quite complex one, since, while the Green's function is the same at both instants (being the reference scenario unchanged), the total fields and are instead different, due to the fact that they implicitly depend on the contrast. However, in the framework we are considering that the variation between and can be assumed to be small (and in case of ischemic strokes also characterized by a weak contrast), so that one can reliably assume that the aforementioned changes in the total fields can be neglected. This corresponds to adopt the Distorted Born approximation and rewrite (2) as follows: where assuming that the adopted reference model is the scenario at the first diagnosis, expresses the variation of the stroke region with respect to the diagnosed one.
The imaging problem we have to face is hence cast as the reconstruction of the spatial map of , from the measured differential scattered field . Notably, thanks to the above considerations, this problem is a linear one, since the relationship between the data and the unknowns is expressed through the linear operator . This is an important circumstance, as the solution of linear ill-posed inverse problems is free from occurrence of false solutions. On the other hand, the operator is compact, owing to the properties of its kernel [10]. As such, the imaging task, in order to provide reliable diagnostic results, has to be carried out adopting suitable regularization tools to counteract the ill-posedness.
In doing so, one has also to consider that a quasi-real-time follow-up monitoring is the ultimate goal, in order to promptly detect sudden changes, so that the adopted imaging procedure should rely on fast and simple methods, able to provide reliable images in a possibly short time. This issue is addressed in the next section.
3. Inversion Procedure: Exploiting an Approximated Reference Scenario
According to the formulation of the monitoring problem derived in the previous Section, the imaging problem at hand is cast through a linear and ill-posed inverse problem. As well known, for this class of problems, a reliable and flexible regularized inversion tool is the TSVD scheme [10]. In such noniterative scheme, the sought differential contrast is explicitly estimated through the inversion formula: where , , and denote the singular values and the left and right singular vectors of the linear operator , respectively. In (4), the truncation index acts as regularization parameter and its value has to meet the trade-off between stability of the reconstruction (with respect to noise on data and modeling errors arising from the adopted reference scenario) and accuracy of the resulting image (in terms of spatial resolution). If the noise level is known, a suitable value could be obtained by using the Morozov discrepancy principle or the L-curve method [10]. However, since in practice such a knowledge is not available, we set the truncation index to a value which is slightly lower than the one corresponding to the change of slope of the singular value's spectrum (theoretically foreseen by the analytic nature of the operator's kernel). It is also worth to note that, in the clinical problem we are considering that this issue could be less important than that in other contexts. As a matter of fact, when performing the examination, the user knows where the variation is expected to occur, so that she/he can discard possible artifacts arising in other areas of the brain (that may result from setting the threshold beyond the suitable value). Thus, one can think that in practice, thanks to the real time nature of the image formation procedure, the clinical operator can adaptively “fine tune” the threshold value with respect to the one provided by the algorithm.
The reference scenario is required in the inversion to build the kernel of the integral operator (through the Green's function and the corresponding total field). In the ideal case, the inner structure of the head at the time of the first diagnosis would be exploited to this end. Unfortunately, an exact knowledge of such a scenario is never available, so that some strategies have to be devised to overcome this drawback.
To this end, a first (optimistic) possibility is to assume that a not perfect, but at least accurate, description of the scenario at the time is available, for instance, via the conversion of a high resolution diagnostic image into relative permittivity and conductivity maps. However, such a translation would be unavoidably affected by uncertainties, as the actual electromagnetic properties of different human tissues change from patient to patient. Moreover, a proper exploitation of a diagnostic image provided by another technique requires to tackle the nontrivial image registration problem, that is, the process of transforming different sets of data into a unique coordinate system.
An alternative, and much more simple solution is suggested by the linear nature of the imaging problem at hand. As a matter of fact, within a linear framework, the nature of the scattering phenomenon (i.e., the field's perturbation due to the contrast variation) is inherently local, so that the imaging performances are expected to be mostly influenced by the properties of the tissues surrounding the stroke, rather than those of the overall scenario (in our case the entire brain). Since the main kind of tissues constituting the brain (in the area where strokes typically occur) is white matter, a possible strategy amounts to build the required kernel using a simple reference model, in which the head is “filled” with a homogeneous medium having the average electric properties of white matter.
Accordingly, in the following numerical analysis, we have appraised the method's performance in both ideal conditions (i.e., when an almost exact map of the brain tissue electric properties is available) and for the simple, nonspecific, scenario described above. It is important to remark that a meaningful difference exists between the two cases from a computational point of view, since using the nonspecific scenario, allows us not only to overcome the aforementioned issues, but also to achieve images in real-time. As a matter of fact, the computationally intensive parts of the TSVD image formation procedure consist in the construction (through the solution of several forward scattering problems) of the matrix representing the discretized version of the operator and the evaluation of its singular system. In principle, this matrix depends on the adopted electromagnetic scenario and therefore should change from patient to patient. Conversely, the possibility of performing this computation with respect to a non-patient-specific scenario allows moving this task to an offline stage, thus only limiting the actual image formation stage to the actual evaluation of (4). This latter is a computationally inexpensive task, as it requires to evaluate scalar products involving vectors of size .
In the following numerical analysis, we will check the validity and the performance achieved when applying the proposed differential TSVD imaging scheme using this simplified scenario to build the inversion kernel.
In Figure 1, we report the flow chart of the proposed procedure, in which the role of the reference scenario is highlighted.
4. Numerical Analysis
To assess the achievable reconstruction capabilities of the proposed strategy, we have considered the 2D reference permittivity profile shown in Figure 2(a), which corresponds to the slice number 63 (in the direction) of the Zubal anthropomorphic head phantom [19]. The dielectric properties of the involved biological tissues have been set according to [20]. An ischemic stroke has been simulated as an elliptic inclusion positioned at cm and cm with major and minor semiaxis of cm and cm, respectively, and having electric properties and S/m.
(a)
(b)
The working conditions have been set according to the guidelines reported in our previous studies [14]. In particular, these results suggest that frequencies lower than 1 GHz represent a suitable choice. Hence, to accommodate the trade-off between wave penetration and achievable spatial resolution, we assume that probes work at 1 GHz and that the head is supposed to be surrounded by a matching medium having electric properties and S/m, which can be for instance achieved through a proper mixture of water and Triton X100 [21]. Note that the use of a matching liquid requires the adoption of suitable solution when setting the device in practice, as done for instance in [18], where antennas are covered with plastic bags filled with the adopted matching medium.
In the simulations, the differential data have been gathered using infinite length filamentary antennas oriented along the -axis. Sources and measurement points are evenly spaced on a circumference having radius cm surrounding the investigated domain, that is, a square of side cm. Note that the number of probes herein considered is chosen in such a way to collect all the available information, according to the Degrees of Freedom Theory (DOF) [22, 23].
The first reported example deals with the followup of an ischemic stroke. In particular, in Figure 2(a), the stroke in its initial stage is shown, that is, at the instant corresponding to the first diagnosis, while the black contours describe its evolution in time, at the instants , , and . In Figure 2(b), it is shown the nonspecific reference model adopted to perform the inversion. For the sake of comparison, also an almost ideal case in which the reference model has been built by corrupting the exact relative permittivity and conductivity maps shown in Figure 2(a) with a ±5% random error on each pixel. The results of this first comparison are shown in Figures 3(a)–3(f). As can be seen, the lack of patient specific information does not significantly affect the obtained images and, more importantly, still allows to follow the evolution of the stroke over time. In the ideal case, the TSVD truncation index has been set equal to , while in the simplified reference scenario it is .
(a)
(b)
(c)
(d)
(e)
(f)
The second example deals with the effect of a mismatch in the position of the probes between two consecutive measurement times. To this end, considering that the measurement apparatus is typically constituted by a fixed geometry array, so that probes can not change their relative position, a displacement of degrees has been introduced. Note that in this case the almost ideal reference scenario has been adopted (). As can be seen from Figure 4, in this situation, the adopted TSVD based approach still provides a reliable reconstruction. In particular, the inaccurate positioning of the device results in some artifacts located in the corners of the imaged domain. While this is of course a not desirable outcome, the fact that these inaccuracies exhibit symmetries entails the possibility of recognizing them as spurious artifacts, possibly suggesting the medical staff the repositioning of the apparatus.
(a)
(b)
(c)
The final assessment has been carried out in almost realistic conditions, that is, using the nonspecific reference scenario and assuming that data are affect by both an erroneous positioning of the probes (an offset of 2 degrees) and measurement noise. In particular, synthetic data have been corrupted by a gaussian noise ( dB with respect to the total field). It is worth to note that such a noise level already entails that the useful signal (the differential scattered field) is almost overwhelmed by noise. This means that one has to take this circumstance into account in designing the device, as a constraint on the system's dynamic range. In this respect, it is certainly interesting that such a noise level is already achieved in MWI systems for other medical applications. Figure 5 shows the TSVD images (for ) and confirm the capability of the proposed approach of correctly tracking the evolution of the stroke, even in the considered realistic conditions. This can be ascribed to the intrinsically low-pass nature of the TSVD procedure, which is indeed robust against the presence of uncertainties on the reference scenario, as well as on the measured data.
(a)
(b)
(c)
In order to quantitatively appraise the above results, we have computed the relative root mean square (RMS) error of the normalized images in a region of interest surrounding the stroke. In particular, we have considered a rectangular area centered on the stroke 12 cm × 8 cm wide. Table 1 summarizes the error values obtained in all cases, which fully confirm the conclusion that the approach is capable to provide accurate results without requiring patient-specific information.
As a further example, we have simulated the followup of a hemorrhagic stroke during the reabsorption stage, which is of interest to monitor the efficacy of drugs therapies. To this end, we have considered the phantom corresponding to slice 45 of the Zubal phantom [19] and the dimensions of the elliptic shaped stroke at times , , , and are, respectively, cm and cm, cm and cm, cm and cm, and cm and cm. In this case, the electric properties of the stroke are those of the blood at GHz (, S/m). Corresponding permittivity profiles are shown in Figure 6(a), the obtained reconstructions (for ) are reported in Figures 6(b)–6(d). Again, the RMS error assumes quite low and stable values (RMS, RMS, and RMS).
(a)
(b)
(c)
(d)
5. A 3D Validation
The previous results were related to a 2D geometry. Clearly, such a schematization is only a simplification of the actual scenario and cannot fully account for the complex (vectorial) interaction between electromagnetic waves and human tissues. Accordingly, in the following, we report a 3D first validation of the above concepts, in order to extend the findings achieved in 2D also in the 3D case.
When facing the 3D problem, the increase of the number of probes required to correctly acquire the information poses an issue concerned with the overall complexity of the apparatus and the required measurement time. For this reason and taking into account that DOF in the 3D case provides only an (loose) upper bound to the actually available information [23], the proof-of-concept reported in the following is aimed to show the feasibility of the proposed follow-up procedure, with respect to an apparatus whose overall number of antennas is in the same order as the “ideal” one adopted in the 2D case. Accordingly, we have considered a simple measurement configuration, in which elementary electric dipoles are displaced on three circles of radius cm. On each circle evenly spaced dipoles are located, but the relative position of the probes between the circles is not the same, but shifted of . As a matter of fact, such a staggering of the antennas allows us to increase the spatial diversity among the different measurements, thus increasing the amount of collectable information [24].
Another issue that has to be tackled in 3D is concerned with the required computational burden. In this respect, the above demonstrated capability of computing the inversion kernel with respect to a simplified nonspecific model is a crucial outcome, as it allows us to still build the image in real-time, performing the computationally intensive part of the algorithm offline, even before the measurement stage. On the other hand, this offline stage becomes obviously more demanding with respect to the 2D case, especially as far as the requirements of the computational facilities in terms of RAM are concerned. For this reason, to make the overall simulation feasible using a standard 8 GB RAM computer, but without affecting in any way the significance of the validation, we have considered a scaled version of the anthropomorphic Zubal phantom [19], in which the head has been reduced to cm (instead of cm), by scaling all the quantities.
The simulated ischemic stroke is supposed to lie at cm, cm, and cm. At time its size is cm, cm, and cm; at time its size is cm, cm, and cm; at time its size is cm, cm, and cm. In Figure 7(a) the slice cutting the center of the stroke along the -axis is shown, with the black contours indicating the evolution of the stroke area. The reference model has been built with respect to the portion of the head enclosed by the probes and the synthetic data have been corrupted by a gaussian noise with SNR = 70 dB (with respect to the total field). The TSVD truncation index is . The , , and sections of the TSVD reconstructions across the stroke at time and are reported in Figures 7(b)–7(d) and Figures 7(c)–7(e), respectively. The RMS error as computed in a volume of size 8 × 4 × 5 cm centered on the stroke is at and at . As can be seen, reliable images are still obtained using a simple measurement apparatus close to actually feasible systems [18].
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
6. Conclusions
We have proposed and assessed a robust procedure for brain stroke followup, which is an application of emerging clinical importance in stroke treatment [3]. In particular, we have shown that an inversion strategy based on a TSVD algorithm is capable of obtaining reliable and readable images of the stroke's evolution even without patient-specific information and in presence of inaccuracies in the measurement process, arising from positioning errors from one acquisition to the other and measurement noise.
In particular, the possibility of achieving the diagnostic image without translating MRI or CT images into electric maps of the brain, but rather exploiting a very simple nonspecific reference model, is crucial to have an actual possibility of practical exploitation. As a matter of fact, the computationally intensive part of the processing scheme is pursued offline, while the actual image formation is carried out in real time.
Finally, the preliminary assessment in the realistic 3D scenario has shown that the procedure can work even with very simple measurement configurations. Accordingly, ongoing efforts are focused on the optimization of the apparatus in order to obtain better images, while keeping the system as simple as possible. For example, an alternative way to exploit a simple measurement apparatus like the considered one is moving the system and collecting the scattered field at different quotas.
Conflict of Interests
The authors declare that there is no conflict of interests.