1. science
2. mathematics

UNIVERSITY OF PAVIA
Department of Electrical Engineering
Dottorato di ricerca in Ingegneria Elettronica, Informatica ed Elettrica
Doctorate dissertation by Maria Evelina Mognaschi, MEng Inverse Problems in Bioelectromagnetism: Experimental and Theoretical Aspects Supervisor Professor Antonio Savini Academic year 2005/2006 To my dad Contents Introduction Chapter I ‐ INVERSE PROBLEMS IN BIOELECTROMAGNETISM
I.1 Introduction I.2 Direct problems I.3 Inverse problems Chapter II ‐ CASE STUDY 1: NON‐INVASIVE THERMOMETRY FOR THE THERMAL ABLATION OF LIVER TUMOR
II.1 Introduction II.2 Fundamentals II.2.1 Problem statement II.2.2 Solution Strategy II.2.3 Case study II.2.4 Numerical results II.3 Localized perfusion model II.3.1 Problem statement II.3.2 Case study II.3.3 Numerical results II.4 Conclusion Chapter III ‐ CASE STUDY 2: SYNTHESIS OF COILS FOR MAGNETIC SIMULATION OF PERIPHERAL NERVES
III.1 Introduction III.2 Experiments III.2.1 The collision method revisited: procedure III.2.2 The collision method revisited: results III.2.3 Measurements of the time constant of a nervous fiber with different transducers: procedure
III.2.4 Measurements of the time constant of a nervous fiber with different transducers: results III.3 Automatic procedure for the synthesis of the optimal coil III.3.1 The direct problem 1
4
5
7
13
16
17
19
19
21
22
26
28
29
31
32
36
37
38
40
42
45
48
49
54
54
III.3.2 The model III.3.3 Inverse problem III.3.4 Numerical results III.4 Conclusion Conclusion Bibliography Papers authored by the candidate 60
66
69
80
81
84
92
Introduction 1
The branch of engineering called bioengineering is a bridge between medicine and engineering. On one side, physicians do not have a suitable physical and mathematical background to model some biomedical problems; on the other side, engineers do not have enough medical background to understand data or sufficient practise to interact with a patient’s body for experiments. For years there has been close collaboration between physicians and engineers and every year this collaboration grows more and more. The need for this interaction is due to the complexity of the human body and in general of every biological mechanism of living beings. This merge of knowledge of biology, medicine and electrical engineering in particular is also reflected by congresses and journals. In fact, most congresses whose topics had been traditionally numerical methods in classical electromagnetism, recently dedicated one or more sessions to bioelectromagnetism. This is the case of OIPE (Optimisation and Inverse Problem in Electromagnetism) which in 2006 had among its topics the biomedical systems. Similarly, ISEM, which in 2005 changed its name from “International Symposium on Applied Electromagnetics and Mechanics” in “International Symposium on Interdisciplinary Electromagnetic, Mechanic and Biomedical Problems”, in 2005 had a session called “Bioelectricity and Biomagnetism”. The last IGTE in Graz (Symposium on Numerical Field Calculation in Electrical Engineering, 2006) proposed biomedical applications in poster and oral sessions. Also large‐audience congresses in electromagnetism are including biomedical topics. An example comes from CEFC 2006 (IEEE Conference on Electromagnetic Field Computation) in which an oral session was on “Nanomagnetics and Biomedical Application” and a poster session on “Bioelectric Field Computation”. Other congresses that targeted on the bioelectromagnetic topics are the IEEE International Conference on Bioelectromagnetism or the IEEE Bioengineering Conference. Moreover, journals on classical electromagnetic problems, such as COMPEL (The International Journal for Computation and Mathematics in Electrical and Electronic Engineering), IJAEM (International Journal of Applied Electromagnetics and Mechanics), IEEE Transactions on Magnetics, are now publishing more and more works on bioelectromagnetic problems. 2
All this necessity to publish or to discuss bioelectromagnetic problems is fed by numerous research groups all over the world. Just to remember some groups of people working on inverse problems in bioelectromagnetism, the research group of Prof. Brauer, The Biomagnetic Group of Ilmenau (http://www4.tu‐
ilmenau.de/EI/ATE/bio/biomag_tui.html) has been working for many years on inverse problems in bioelectromagnetism; in Finland there is the Finnish Centre of Excellence in Inverse Problem Research (http://math.tkk.fi/inverse‐
coe/index.html) and also the Ragnar Granit Institute of the University of Technology of Tampere, the director of which is Prof. Malmivuo (http://www.rgi.tut.fi/www/index.htm). In my research group, in the CAD laboratory of the Electrical Engineering Department of the University of Pavia, inverse problems in quasi‐static electromagnetics have been traditionally studied, but now attention is also given to biomedical applications. My PhD thesis comes out from this rapidly changing environment and it is set in this framework of interest. In Chapter I an introduction to direct and inverse problems in bioelectromagnetism is given. The following two chapters are related to the two case studies considered, respectively. In Chapter II a non‐invasive thermometry method to be applied during a treatment of thermal ablation of liver cancer is presented. This case study was in collaboration with the Department of Electrical Engineering of Padua. In Chapter III a synthesis of a field source for electromagnetic stimulation of peripheral nerves is performed. In Chapter III some experiments, in collaboration with the “Istituto Neurologico C. Mondino” of Pavia, are also presented. The collaboration with Prof. P. Di Barba for all studies presented in this work is gratefully acknowledged. 3
Chapter I INVERSE PROBLEMS IN BIOELECTROMAGNETISM
4
I.1 Introduction The aim of collecting data is to gain meaningful information about a physical system or phenomenon of interest. However, in many situations the quantities to be determined are different from the ones which are possible to be measured. If the measured data depends, in some way, on the quantities to determine, then the data at least contains some information about those quantities. Starting with the data measured, the problem of trying to reconstruct the quantities wanted is called an inverse problem. Loosely speaking, an inverse problem is one in which an effect is measured and the cause of it is to be determined. Inverse problems for which data come from measurements are known as identification problems. Another group of inverse problems refers to synthesis problems in which data are assumed arbitrarily. To this group belong optimal design problems, where the purpose is to design a device which can provide a defined behaviour or optimal performance (Neittaanmäki et al., 1996). The first group of inverse problems is very popular in bioelectromagnetism. All the problems in which a diagnostic image has to be reconstructed belong to this category. For instance, the computer axial tomography is a technique to obtain slice images through the body in vivo. It is known that X‐rays are partially transmitted through the body, and that the opacity of various internal structures to X‐rays varies, so that a picture of the variation of the absorption coefficient in the body would give a good picture. However, the only measurements that one can make non‐invasively is to shine X‐rays through the patient and to measure the total absorption along lines through the body. Given a collection of such line integrals (the “data”), the problem is to reconstruct the absorption as a function of position in the body (the “image”). To this group of inverse problems also belong the identification of electromagnetic sources inside the human body. In fact, some tissues and organs of the human body are able to produce electromagnetic signals; above all the brain and the heart. From here comes the possibility to reconstruct the sources of these signals, often recorded on the body surface. Typical examples 5
are the electro‐ and magnetoencephalography (Küçükaltun‐Yıldırım et al., 2006) and electro‐ and magnetocardiography techniques (Brauer et al., 2000). The first case study reported in this work belongs to this kind of inverse problems. In fact, the problem is to reconstruct, in a non‐invasive way, the temperature field on a slice of a patient’s torso when the tissues are heated from an external source and the known data are the temperatures recorded on the patient’s skin. In general, solving these problems can give us knowledge of inner parameters, geometries, tissue properties and field distributions inside human beings, without invasive measurements. In turn, optimal design problems are very popular in classical engineering. In the biomedical field they are applied to the design of electromedical devices (from MRI machines (Cavaliere et al., 2000) to coronal stunts), transducers and prostheses (Ngale Haulin and Vinet, 2003). The second case study presented in this work belongs to this category of inverse problems. A device for magnetic stimulation of nervous tissues is taken into account. In particular the shape of the coil, in which a time‐varying current flows, is modified in order to have a selective and powerful stimulation of peripheral nerves. The solution of some kinds of inverse problems, e.g. optimal design problems, requires a module for calculating the field. This module is associated with an optimization module which performs the minimization of some suitable objective functions. In traditional computer‐aided design these two modules are interconnected in a way which simulates a trial and error procedure. So, starting from an initial design, the field is analysed by means of analytical, semi‐analytical or numerical methods. Then, a check is made as to whether the device has the desired properties or the objective functions are minimized sufficiently, according to a prescribed tolerance. If it does not, some parameters of the problem are varied, the field is recalculated again and so on. At each iteration at least one field solution has to be calculated. This procedure is quite troublesome and time consuming; a compromise between accuracy and consuming time has to be made. 6
I.2 Direct problems From an engineering point of view, the starting point in studying biomedical problems is to implement a model of the system to be studied. For these kinds of problems, this is not a simple task. In fact, the properties of the object to be modelled are subject to great variability between individuals and between different conditions of the subject, such as body or external temperature, anxiety state, hyper or hypo‐ventilation state and so on. Moreover, in bioelectromagnetism the geometry of the direct problem can be even very complicated. If the problem is to be solved numerically, the domain has to be built. For simple geometry, the domain can be drawn by hand with the software used for the numerical solution. This is the case of certain human body parts, such as arms, legs (but not knees or elbows) or subparts of organs such as the windpipe, the oesophagus or blood vessels. Entire organs e.g. the brain or slice of torso present a complicated geometry. Because of the difficulty to obtain a good accuracy in drawing this kind of domains by hand, some specific softwares have been developed (Ziółkowski and Brauer., 1996; www.amiravis.com). These softwares allow the users to look at the medical images e.g. in Dicom format and to modify these images in order to highlight the anatomical areas of interest and to cancel the non‐interesting parts. These anatomical parts could be small fat isles, little bone excrescences or imaging artefacts; the choice of the parts to eliminate depends, obviously, on the purpose of the solution. Moreover, most of these softwares are able to generate the discretization of the domain and to pass these objects (domain and grid) to other commercial softwares for numerical solution. The inverse problems proposed in this work are based on the numerical solution of the direct problem. In particular the finite element method is used. The direct problem of the first case study consists of field calculations on a slice of a patient’s torso. For the reasons explained above, a CT image is considered (see Fig. I.1) and elaborated with the commercial software Amira. 7
Backbone
Liver Muscle
Fig. I.1 – CT image of a patience torso from which the domain of the first case study is imported. Then, the domain is passed to a finite element simulator and, after assigning the material properties, the boundary conditions and the sources, it is solved. The practical operations concerning the direct problem of this case study were implemented by researchers at the Department of Electrical Engineering of Padua. A different strategy for the second case study is used. It was necessary to model an arm in three dimensions. Because of the cylindrical symmetry of the part and the minor number of different tissues it was simpler to draw it by hand than to import it by some pre‐processor software. In Fig. I.2 a MRI image of a forearm is represented. 8
Muscle
Ulnar nerve Bone
Fat Fig. I.2 – MRI image of a forearm: the geometry is simple. After drawing the forearm, the domain so built was passed to a finite element simulator. The problem to be solved is a time‐varying problem in three dimensions. From the mathematical point of view it belongs to the diffusion problems. A method used for solving 3D finite element problems is a combined electric vector potential T and scalar magnetic potential Ω approach, commonly known as the T ‐Ω method. The T ‐Ω method for diffusion problems is here described. Starting from the Maxwell equations for the quasi‐stationary approximation in a conducting medium: ∇⋅D = ρ (I.1) ∂B
∇×E = −
(I.2) ∂t
9
∇⋅B = 0 (I.3) ∇×H = J (I.4) where D is the dielectric induction, E the electric field, B the magnetic induction, H the magnetic field, J = J 0 + σE with J 0 the impressed current density, σ the electric conductivity and ρ the volume charge density, it is possible to introduce the electric vector potential ( T ) and the scalar magnetic potential (Ω).
In fact, from Eq. (I.4) and knowing that the curl of a scalar quantity is null, it comes that J = ∇ × H = ∇ × ( H + ∇ Ω) = ∇ × T (I.5) from which (I.6) H = T − ∇Ω In order to define T uniquely, a gauge for T has to be introduced. The diffusion equation, governing the quasi‐static magnetic problems, can be expressed in terms of the two potentials: ∂
∇ × (σ −1 ∇ × T ) = ∇ × σ −1 J 0 − µ( T − ∇ Ω) (I.7) ∂t
where µ is the magnetic permeability. Assuming J 0 =0 and imposing the gauge ∇⋅T =
∂Ω
µσ ∂t
(I.8) two independent equations for T and Ω are obtained: 10
∇ 2 T − µσ
∂T
=0 ∂t
(I.9) ∇ 2 Ω − µσ
∂Ω
= 0 ∂t
(I.10) Solving Eq. (I.7) is equivalent to minimizing an energy functional derived from the Poynting theorem. Taking into account Eq. (I.2) and a fundamental vector identity: ∂B
∇ ⋅ ( E × H) = H ⋅ (∇ × E) − E ⋅ (∇ × H) = − H
−E⋅J (I.11) ∂t
Referring to a simply‐connected domain D bounded by a surface S in the time interval [0, T], the Poynting theorem holds: 1 ∂
(µH ⋅ H) dD − ∫ E ⋅ J dD (I.12) ∫ (E × H) ⋅ n dS = ∫ −
2 ∂t
S
D
D
Substituting Eq. (I.6) and the constitutive equation J = σE in Eq. (I.12) gives: 1 ∂
[µ(T − ∇Ω) ⋅ (T − ∇Ω)] dD +
2 ∂t
D
−1
− ∫ σ ( J ⋅ J) dD
−1
∫ σ [ J × (T − ∇Ω)] ⋅ n dS = ∫ −
S
(I.13) D
Integrating over time, one has: T
1
−1
∫ ∫ σ [ J × (T − ∇Ω)] ⋅ n dt dS = ∫ − [µ(T − ∇Ω) ⋅ (T − ∇Ω)] dD +
2
S0
D
T
(I.14) −1
− ∫ ∫ σ ( J ⋅ J) dt dD
D0
11
From Eq. (I.14), the energy functional is derived: T
1
χ( T , Ω) = ∫ ∫ σ −1 J × ( T − ∇Ω) ⋅ n dt dS + ∫ µ ( T − ∇Ω) ⋅ ( T − ∇Ω) dD
2D
S0
(I.15) T
−1
+ ∫ ∫ σ ( J ⋅ J ) dtdD
[
]
[
]
D0
The first term is related to the energy delivered out of region D, the second to the stored magnetic energy and the third to the Joule losses (Preston and Reece, 1982). Setting the variation of χ(T − ∇Ω) , with respect to T and to Ω, to zero, Eqs. (I.9) and (I.10) are found, respectively. Hence, solving Eqs. (I.9) and (I.10) is equivalent to minimizing the functional of Eq. (I.15), with T and Ω regular enough to have all the terms of Eq. (I.15) non‐
singular. Once T and Ω are calculated, it is possible to obtain the electric field E . In fact, referring to Eq. (I.5), one has: J = J 0 + σE = ∇ × T (I.16) Hence E = σ −1 ( ∇ × T − J 0 ) (I.17) This method is not able to treat a multiple connecting conductor. If this were the case, some strategies could be used e.g. filling the hole with a low‐
conducting material. However, the domain used for the second case study is simply‐connected. The great advantage of the T ‐Ω method is that in non‐conducting regions, it is necessary to solve only a problem governed by Eq. (I.10). In such a way, only a variable is associated to each node and the number of degrees of freedom are 12
decreased with respect to the full T ‐Ω problem (Eq. (I.7)). The software used for solving the finite element problem of the second case study proposed is MagNet, a commercial software by Infolytica, a Canadian software house (www.infolytica.com). For 3D problems, MagNet implements the T ‐Ω method. I.3 Inverse problems
Inverse problems are usually solved by means of a minimization of a suitable error function. Sometimes the choice of the error function is strategic and very important. In fact, the data measured or simulated have to be related to the unknown variables of the problem, preferably in a direct and close way. If there is a lack of information or the error functional is badly choosen the inverse problem can be ill‐posed in the Hadamard sense (Hadamard, 1923), i.e.: 1. Existence of the solution. 2. Uniqueness of the solution. 3. Stability of the solution, i.e. the solution depends continuously upon the data. More precisely, if the problem is determined by an equation Kf=g, K:X→Y with X and Y normed spaces, then the problem is well posed if K is bijective and if the inverse operator K‐1:Y→X is continuous (Neittaanmäki et al., 1996). In bioelectromagnetism, inverse problems are often ill‐posed. Let’s think, for instance, of the electrocardiography. Its aim is to restore the heart activity from a given set of body surface potential maps. It is an ill‐posed inverse problem. In fact, the solution is unstable i.e. an arbitrary small change of the source signal can cause an arbitrary large change of the solution. The aim of this work is not to study the inverse problems presented from a theoretical point of view; however some considerations can be made. For the inverse problems described in the further chapters, the existence of solutions is expected. In fact, for the physics of the problems, their existence is reasonable. 13
Instead, the uniqueness of the solution is not guaranteed. In the case of more than one solution, the use of deterministic methods is inadequate. This is due to the presence of more than one minimum, each one characterized by a different degree of stability. On the other hand, stochastic optimizers are more effective in finding the global optimum because they are independent to the gradient. In this respect, the evolutionary algorithm substantially differs from a deterministic one, in which the search region would be narrowed around the best point in order to converge towards the corresponding, nearest minimum. The drawback, in fact, is that this minimum might be a local one. On the contrary, the evolutionary algorithm, if successful in finding a better point, covers a larger region of search in order to see if there would be another good candidate in the neighborhood. It then does the opposite when this is not deemed possible. This way, there is a non‐zero probability of finding the region where the global optimum of the objective function is located. Moreover, for synthesis problems, the non‐uniqueness of the solution can be an advantage. In fact, in this case the designer can choose between different solutions, which are, in general, different, even if they have the same performances. The stability of the solution is not guaranteed and not expected. In fact it is known that this condition is the most difficult to fulfill in all kinds of inverse problems and in particular in bioelectromagnetism. In the case studies presented, the optimization of the objective function or functions is guided by an evolutionary algorithm of lowest order. For the first case study a single‐objective optimization algorithm is used. In the algorithm used (ESTRA) (Neittaanmäki, 1996) the selection between parent m and offspring x is ruled by the following criterion: the objective function f(x) is evaluated and the test f(x) < f(m) is performed; if the test is successful, m is replaced by x, otherwise m is retained and another offspring is generated. The second case study is characterized by two objective functions; a multi‐
objective optimization algorithm (MOESTRA) is used (Di Barba, 2001). In MOESTRA the selection is ruled by the Pareto criterion, i.e. x replaces m if and only if fi(x) ≤ fi(m) ∀i and fi(x) < fi(m) for at least a value of index i = 1,2. 14
In both the algorithms, the search in the design space begins in a region of radius d0 (standard deviation) centered at the initial point m0 (mean value); m0 is externally provided, while d0 is internally calculated on the basis of the bounds boxing the variation of the design variables. Setting m=m0 and d=d0, the generation of the design vector x = m + u d then proceeds, resorting to a normal sample u ∈ (0 ,1) . It is verified that x fulfils bounds and constraints (i.e. that x is feasible), otherwise a new design vector is generated until it falls inside the feasible region. The selection between parent m and offspring x is ruled by the beforementioned criteria. The next step is concerned with the size of the search region that will be used for the successive iteration. The underlying rationale is that when a point better than the current one is found, the radius of the search region is increased around the new point to search for further improvements; if no improvement is found, the radius of the search region is gradually decreased up to convergence (annealing process). The annealing process is ruled by the history of the k previous iterations, used to establish a trend: if at least a fraction p of the last k iterations were successful (an iteration is successful if x is feasible and improves the objective functions), then the trend is said to be positive, while it is negative otherwise. If the trend is positive, the radius d of the search region is set to q‐1d and otherwise it is set to q d . In particular, during the first k iterations, d remains unchanged. The d
is achieved. Quantities p procedure stops when the prescribed accuracy d0
and q are named probability of success and rate of annealing, respectively and represent the “tuning knobs” of the algorithm; heuristic values for k, p and q are 50, 0.2 and 0.8‐0.9, respectively. The non‐deterministic approach could make the computational cost too high because of the high number of objective function calls, which, in turn, could require the solution of a direct problem. However, the methods used in this work were tested and compared with other algorithms (Di Barba and Mognaschi, 2004; Carcangiu et al., 2006). 15
Chapter II CASE STUDY 1: NON‐INVASIVE THERMOMETRY FOR THE THERMAL ABLATION OF LIVER TUMOR
16
II.1. Introduction The vasculature of a malignant tumor is different from the vasculature of healthy tissue in several ways. It is highly disorganized, angiogenic, leaky and displays extensive branching and shunts. As a consequence, significant regions of a tumor will be nutrition deprived and poorly oxygenated. This means that many solid tumors contain significant fractions of hypoxic and acid regions. Despite the adverse conditions in these regions, they contain viable cells. Such cells are known to be sensitive to the cytotoxic effect of heat, suggesting that heating could be adjuvant in anti‐cancer therapy (Hokland et al., 2006). The procedure of applying moderate, in the sense of nonlethal for cells, thermal doses to a target region is called hyperthermia. The application of lethal thermal doses, resulting in wide spread coagulative necrosis in the affected area, is called thermal ablation. Hyperthermia in oncology has, until now, primarily been used as an adjuvant treatment to conventional therapies, i.e. radio‐ and chemiotherapy, due to the potentiating effect it has to these treatments. Thermal ablation is used as a treatment of primary and metastatic tumor in different parts of the body such as liver, for which the clinical experience is becoming important, lung, kidney and also bone. It has not yet been successfully used for tumors in organs merged in non‐conducting regions such as breast. Energy can be delivered to the target zone in various ways (i.e. radiofrequency, microwave, lasers and ultrasound) and using different devices (i.e. electrodes and antennas). However, up to now, the most used method is the radiofrequency treatment by means of an electrode, implanted percutaneously or laparoscopically in the lesions, under Computer Tomography (CT) or Magnetic Resonance (MR) guidance. The electrode can be needled or hooked, with 2‐5 cm of exposed metallic tip, connected to a radiofrequency generator in the range of 450‐550 kHz. The current flows from the needled electrode to a grounded electrode and the heating effect is mostly due to the electric conductivity of the tissues. From the clinical viewpoint, knowing the internal thermal field is of utmost interest in order to control power deposition resulting in an effective treatment. In fact, the effectiveness consists in reaching the therapeutic temperature (up to 17
42 °C) in the tumoral tissue, preserving, as much as possible, the healthy tissues (0.5 to 1 cm of adjacent tissues). However, the temperature in the target zone must not exceed 105 °C, because in this case there would be bubble formation, boiling of tissue and eventually carbonization (Ulucakli, 2006). A carbonized tissue presents a high resistivity; in this way the heating process could be blocked. Controlling the temperature is the delicate problem of these kinds of methods. In fact, despite the successes obtained in most hospitals and medical centers all over the world, hyperthermia and thermal ablation treatment have not gained wide‐spread clinical acceptance. This is manly due to the technical inability to control the heating process inside the body of the patient and, in particular, to observe and reconstruct, in a non‐invasive way, the temperature distribution along the target zone and the nearest tissues. In the last years, however, a method to obtain temperature distribution has been improved. It is based on the Magnetic Resonance Imaging (MRI) and consists of acquiring phase images, called baseline images, before heating and subtracting them from phase images acquired during heating. Using standard clinical MRI‐systems typical parameters seem to allow a standard deviation of temperature less than 1 °C with a temporal resolution of about 1 s and a voxel size of 1 mm × 1 mm × 5 mm (Hokland et al., 2006). However, motion artifacts may greatly reduce the accuracy of MRI temperature maps, so respiratory or accidental movement must be corrected (Rieke et al., 2004). Moreover, MRI acquisition is expansive and, above all MRI, cannot be used for individuals with pace‐makers, metallic tips or with prostheses containing electronic devices. Hence, as a possible alternative, it is proposed to identify the temperature distribution in the target volume in a non‐invasive way by acquiring temperatures in a finite number of points on the body surface. To this end, it is necessary to provide a numerical model of the heating process in order to perform the computation of the thermal field, after resolving the electromagnetic problem. A finite‐element model, based on the CT image of a patient, has been developed accordingly. From a mathematical point of view, the problem belongs to the category of continuation problems of the thermal field, governed by the heat conduction equation, starting from experimental data. 18
II.2 Fundamentals The first case considered is based on a multiply connected domain in which the tumor is neglected. The effect of capillaries, uniformly distributed all over the domain, is taken into account. In order to reconstruct the temperature in all the domain, an automatic procedure has been set up. It is based on the thermal finite element (FE) model, of a CT slice of torso. A suitable functional has to be minimized and, at each iteration of the procedure, a thermal FE simulation is performed (D’Ambrosio et al., 2005). To evaluate the functional, the temperature along the thermal probes on the skin has to be known. Hence, preliminarly, a thermal FE simulation, based on the simply connected domain (taking the tumor into account) is solved. The solution so obtained gives the so called “recorded”, “measured” or “estimated” temperature. The solution obtained from the thermal FE solution at each iteration of the identification procedure is the so called “reconstructed” temperature. II.2.1 Problem statement Two 2D domains, Ω and Ω0 are considered (see Fig. II.1). Γ0
Ω0
Ω
Γ1
Γ2
Γ3
Fig. II.1. The domain of the problem. Γ0 ∪ Γ1 ∪ Γ3 is boundary for Ω, while Γ0 is boundary for Ω0. It is supposed that Ω 0 ∩ Ω = Γ0 . Γ0 and Γ1 are continuous boundaries, Γ2 is a discrete boundary 19
consisting of nv points, highlighted in Fig. II.1. In turn, Γ3 is a continuous boundary such that Γ1∩Γ2=0, Γ1∩Γ3=0 and Γ2∩Γ3=Γ2. Along Γ1 the temperature is known; along Γ3 the heat is exchanged by means of convection with the surrounding environment, which is supposed to be at temperature u0.
In view of an application, Γ1 is the electrically grounded electrode, Ω0 is the target volume and Γ2 represents the position of the temperature probes. In a thermal equilibrium state, the direct problem is governed by the bioheat equation (Pennes, 1948) − ∇ ⋅ ( k∇u) = f in Ω ∪ Ω 0 (II.1) where u is the temperature (K), k is the thermal conductivity (Wm‐1K‐1) and f = σE 2 − c b w b ( u − u b ) is the source term (Wm‐3). In turn, σ is the electric conductivity (Sm‐1), E is the known electric field (Vm‐1), cb is the specific heat of blood (J kg‐1K‐1), wb is the mass flow rate (kg s‐1m‐3) and ub is the temperature of blood (K). The term σE2 is the power density due to the electric field, while the term cbwb(u‐ub) represents the cooling effects of the blood flowing in the capillaries (perfusion). Eq. (II.1) is subject to the following boundary conditions: (II.2) u = U along Γ1 ∂u
along Γ3 (II.3) −k
= h (u − u 0 ) ∂n
where h is the convection coefficient and n is the normal unit vector. The direct problem (II.1)‐(II.3), that is assumed to be linear, is a well posed problem; its solution is taken as the estimated temperature field all over the domain Ω ∪ Ω 0 . Conversely, the continuation problem can be formulated as follows: knowing the geometry of the domain Ω, the tissue properties (k, σ, cb, wb, h), the source term f, the boundary conditions along Γ1 and Γ3 and the supplementary condition u * = {u i } i = 1,.., nv along Γ2, find the temperature field u in the Γ2
domain Ω. The solution to the continuation problem is taken as the 20
reconstructed temperature field; ideally, estimated and reconstructed fields should be coincident over Ω; practically, their discrepancy will quantify the accuracy of the reconstructed field. Moreover, the uncertainty of tissue data makes the continuation problem ill‐posed, because it is impossible to guarantee the uniqueness of a solution beforehand (Pałka, 1983; Paruch and Majchrzak, 2004). A remark on the direct problem is worth noting: the process of thermal ablation is dynamic and causes blood coagulation, which reduces perfusion. The perfusion model adopted takes into account the dependence on the temperature and neglects the phase transition for the blood. II.2.2. Solution strategy To solve the problem of continuation numerically, a minimization strategy has been adopted: defining the error functional F(λ ) = u (λ ) Γ − u *
( )
2
Γ2 2
, starting from a guess solution, find λ∗ such that F λ* = inf F(λ ) where u (λ ) Γ is the reconstructed field along Γ2 and u
*
Γ2
λ
2
is the supplementary condition. Hence, the error functional is defined as the two‐norm of the deviation between computed and measured temperatures. In general, λ ∈ R m is a vector of m parameters governing the distribution of the temperature field along Γ0. For the sake of simplicity, λ ∈ ℜ1 has been used, so assuming that Γ0 is isothermal. This assumption is very rough. It has been made just to check if the solution strategy works. However, in the localized perfusion model (Chapter II.3) the isothermal assumption along Γ0 will be removed. Given a value of λ, the direct problem is well‐posed; therefore a finite‐element analysis of the following problem − ∇ ⋅ ( k∇u) = f in Ω (II.4) (II.5) u = U alongΓ1 21
−k
∂u
= h( u − u 0 ) ∂n
along Γ3 (II.6) u = λ
along Γ0 (II.7) is performed and the error functional F(λ) is updated; a derivative‐free minimizer drives the overall procedure in an automated way. The solution strategy is derived from (Colli Franzone, 1980), where a method for the continuation of the electric potential on the epicardial surface was developed. II.2.3 Case study The described solution strategy was applied to a practical situation, during a thermal ablation treatment of liver metastasis (Dughiero and D’Ambrosio, 2004). The 2D domain is a transversal section of the torso passing through the middle of the tumor, directly derived from a CT image of the patient. The domain has been divided into regions representing muscle, fat, liver, tumor and bone. The mesh used to model the domain Ω consists of 2,990 nodes and 5,848 first order‐triangular elements (Fig. II.2). The temperature of the nodes belonging to the grounded electrode is set to U = 37°C; the number of nodes of Γ2 is set equal to the number nv of thermocouples, uniformly placed on the skin of the patient, i.e. nv=10. 22
4
5
Ω0 3
Γ2
6
7
2
8
1
Γ1
9
10
Fig. II.2. Finite‐element model. The needled electrode is placed within the tumoral mass Ω0; the electric potential of the needle is set to 50 V with respect to the grounded electrode Γ1, while the condition of tangential electric field is assumed along Γ3. By knowing the electrical properties of the various tissues, it is possible to solve the Laplace’s equation numerically, and estimate the electric field E and so determine the Joule losses represented by the source term σE2 of the subsequent thermal problem. The electric problem to solve reads: ∇⋅J = 0 (II.8) n⋅J = 0 along Γ3 (II.9) n×J = 0 along Γ1 (II.10) where J = ( σ + jωε 0 ε r )E in the frequency domain. For this problem ωε<<σ, hence, J = σE (Dughiero and D’Ambrosio, 2004). 23
In Fig. II.3 lines of electric potential are represented. Fig. II.3. Electric potential lines. In order to resume the solution procedure, flowchart is reported in Fig. II.4. 24
START
Estimated temperatures (reference solution) Electrical problem on Ω ∪ Ω 0
[Eqs. (II.8), (II.9) and (II.10)]
Thermal problem on Ω ∪ Ω 0
[Eqs. (II.1), (II.2) and (II.3)] Reconstructed temperatures Update λ
Thermal problem on Ω [Eqs. (II.4), (II.5), (II.6) and (II.7)]
Supplementary condition Evaluation of the error function F(λ)
Is F(λ) minimized ?
No
Yes
STOP
Fig. II.4 Flowchart of the numerical procedure. 25
II.2.4 Numerical results After convergence of the minimization procedure, starting from u0=10°C, the following temperature distribution is found: Fig. II.5 – Estimated (pink dot) and reconstructed (black cross) temperatures over the domain. Over a domain D, the global error ε(D) between estimated temperature field uε(i) and reconstructed temperature field ur(i) i = 1,..,np , with np number of nodes of the grid discretizing D, was evaluated in ℓ2 norm: 1
ε( D ) = [( u r − u e ) T ( u r − u e )] 2 (II.11) 26
The result is ε1 ≡ ε(Ω) = 6.121 °C and ε2 ≡ ε( Ω \ Γ0 ) = 1.164 °C. It follows that the error is mainly concentrated along Γ0, due to the assumption of considering Γ0 as an isothermal line. Nevertheless, the constant value obtained for the estimated temperature of Γ0 is a lower approximation to the mean temperature along Γ0 itself. Moreover, because the temperature in Ω0 cannot be lower than the temperature in Γ0, approximating the latter enables one to decide if the given power deposition determines the therapeutic temperature within the tumoral mass. In Table II.1 the temperature recorded by the thermocouples is reported. Table II.1 – Recorded temperature Thermocouple
1
2
3
4
5
6
7
8
9
10
Recorded
temperatures
32.3 31.7 31.4 31.3 31.3 31.5 31.3 31.2 31.3 31.1
(°C)
In Table II.2 the estimated temperatures in n0 = 10 nodes discretizing Γ0 and the reconstructed temperature in the same nodes are reported. Table II.2 – Estimated and reconstructed temperature in Γ0. Estimated Γ0
temperatures 40.2 42.6 44.9 42.2 43 40.1 41.1 42.3 42.6 42.7
(°C)
Reconstructed
40.9
Γ0 temperature
(°C)
In Table II.3 global errors and reconstructed temperature of Γ0 are reported by varying convection coefficient h and perfusion coefficients w = (w1, w2, w3, w4, w5). It is assumed that h = 12 W m‐2 K‐1, w1 = 16.67 kg s‐1 m‐3 for the liver, w2 = 0.36 kg s‐1 m‐3 for the fat, w3 = 0.45 kg s‐1 m‐3 for the muscles, w4 = 0.416 kg s‐1 m‐3 for the tumor and w5 = 0.54 kg s‐1 m‐3 for the bone (basic problem) (Lang et al., 1999). 27
Table II.3. – Numerical results varying some thermophysical properties. Reconstructed
Variation
ε1 (°C) ε2 (°C)
temperature of Γ0 (°C)
Basic problem
6.121
1.164
40.9
-2
-1
h = 20 (W m °C )
6.132
1.166
40.9
-2
-1
6.142
1.168
40.9
h = 30 (W m °C )
w’ = 0.5 w (kg s-1 m-3) 9.163
2.503
43.9
-1 -3
w’’ = 2 w (kg s m ) 3.662
0.446
39
It can be noted that the reconstructed temperature is somewhat sensitive to the considered variations. In Table II.4 global errors, reconstructed temperature of Γ0
and the residual of the error functional are shown for different starting points of the minimization procedure. Table II.4. – Depedence of numerical results on the starting point of the minimization. Starting
Reconstructed Final residual
of the error
temperature of ε1 (°C) ε2 (°C) temperature of
functional (°C)
Γ0 (°C)
Γ0 (°C)
10
6.121
1.164
40.8951
1.5827 10-4
30
6.107
1.161
40.8952
1.7150 10-4
70
6.107
1.161
40.8952
1.7118 10-4
1.5831 10-4
6.12
1.164
40.8951
95
It can be noted that the solution to the continuation problem is stable with respect to the initial guess for a given accuracy. II.3 Localized perfusion model Starting from the basic problem, an improved direct problem is considered. The perfusion due to big vessels is now considered, as well as the perfusion due to capillaries (D’Ambrosio et al., 2006). Whenever large vessels are close to the tumor, their effect on heat exchange has to be considered (Tungjitkusolmun et 28
al., 2002; Ulucakli, 2006). In fact, the blood flowing in the vessels takes away heat from the region under heating, decreasing the thermal effect. II.3.1 Problem statement In this case three 2D domains Ω, Ω0 and Ω1 are considered (see Fig. II.6). Γ0
Ω0
Ω1
Ω
y x Γ1
Γ2
Γ3
Γ4
Fig. II.6 – Field domain: cross section of a human abdomen. Γ0, Γ1, Γ2 and Γ3 are continuous boundaries, while Γ4 is a discrete set of n points, highlighted in Fig. II.6, such that Γ3∩Γ4=Γ4. It is supposed that Ω ∩ Ω 0 = Γ0 and Ω ∩ Ω 1 = Γ1 . In view of the simulation of the thermal ablation treatment, Γ2 is the grounded electrode, Ω0 is the target zone, Ω1 is a blood vessel, Γ3 is the skin surface and Γ4 represents the position of the temperature probes. In a thermal equilibrium state, the direct problem is governed in {Ω ∪ Ω0 } by the bioheat equation: − ∇ ⋅ ( k∇u( x , y)) = σE 2 ( x , y) − c b w b ( u − u b ) (II.12) subject to the following boundary conditions: u = U ∂u
−k
= h 1 (u − u 0 ) ∂n
along Γ2 (II.13) along Γ1 (II.14) 29
−k
∂u
= h 2 (u − u b ) ∂n
along Γ3 (II.15) where h1 is the convection coefficient for the skin, while h2 is the convection coefficient for the vessel wall, u0 and ub are the temperatures of air and blood, respectively. The inverse problem then reads : knowing the geometry of the domain Ω, the tissue properties (k, σ, cb, wb, h1, h2), the source term f, the boundary conditions along Γ1, Γ2 and Γ3 and the supplementary condition u * Γ = {u i } i = 1,.., nv 4
along Γ4, find the temperature field u in the domain Ω. The error functional to be minimized is F(λ ) = u(λ ) Γ4 − u *
At each iteration, the direct problem to be solved is: − ∇ ⋅ ( k∇u) = f u = U ∂u
−k
= h 1 ( u − u 0 ) ∂n
∂u
−k
= h 2 ( u − u 0 ) ∂n
u = λ
Γ4 2
, with λ ∈ R m . in Ω (II.16) alongΓ2 (II.17) along Γ1 (II.18) along Γ3 (II.19) along Γ0 (II.20) 30
II.3.2 Case study As for the previous problem, a CT slice of a torso has been discretized with triangular elements. The domain has been divided into regions representing muscle, fat, liver, tumor, vessels and bone. The mesh used to model the domain Ω consists of 11,935 nodes and 23,608 first order‐triangular elements (Fig. II.7). Fig. II.7 – Finite element model. For the grounded electrode U = 37°C, for the convection exchanges h1=800 Wm-2K-1 and ub=37 °C, h2=12 Wm-2K-1 and u0=25 °C; the number of nodes of Γ4 is set equal to the number nv of thermocouples, uniformly placed on the skin of the patient, i.e. nv=26. The electric problem to solve is described in Eq.(II.8), subject to the following boundary conditions: n⋅J = 0 along Γ3 (II.21) n× J = 0 along Γ2 (II.22) 31
In Fig. II.8 lines of electric potential are represented. Fig. II.8 – Electric potential lines II.3.3 Numerical results
First, a constant temperature along Γ0 is considered ( λ ∈ ℜ1 ). The starting point is 50°C and the range of variation for λ is [0;100] °C. The final point of the minimization is λ*=46.1 °C while the mean of the temperature estimated along Γ0 is 48.7 °C. Graphical results are shown in Fig. II.9. 32
Fig. II.9 – Estimated (pink dot) and reconstructed (black cross) temperature over the domain. The region around the tumour and the vessels is the most critical. In fact, in that region, there is a large temperature gradient; this contrasts the assumption of constant temperature along Γ0. In Fig. II.10 a detail of the zone is represented: 33
Fig. II.10 – Estimated (pink circle) and reconstructed (black cross) temperature over the boundaries of tumor and vessels. The errors ε1 ≡ ε(Ω) and ε2 ≡ ε( Ω \ Γ0 ), with ε defined as in Eq.(II.11), are calculated: ε1 = 21.581 °C and ε2 = 12.291 °C. In the region closest to the tumor, discrepancies between estimated and reconstructed temperature are found, however the solution is acceptable, because the reconstructed temperature is near to the mean of the estimated temperature along Γ0. Moreover, the temperatures found under‐estimate the mean value along Γ0 and this better fits our purpose. In order to minimize the errors ε1 and ε2, a two‐variable linear interpolation of the temperatures along Γ0 is introduced: λ ∈ ℜ 3 , i.e. λ=[ λ1, λ2, λ3], for which u r Γ = λ 1 x + λ 2 y + λ 3 , is now considered. 0
The following limits for λ are considered: ‐500<λ1<0, ‐500<λ2<0 and 0<λ3<500. The starting point is λ=[0 0 40]. The minimization procedure gives λ*=[37.5 ‐374 101] and ε1 = 13 °C and ε2 = 7.7 °C. The results are shown in Fig. II.11 and in Fig. II.12. 34
Fig. II.11 – Estimated (pink dot) and reconstructed (black cross) temperature over the domain. Fig. II.12 – Estimated (pink circle) and reconstructed (black cross) temperature over the boundaries of tumor and vessels. 35
Both of the errors decrease with respect to the constant‐interpolation case. The mean of the reconstructed temperatures is 49.4 °C, closest to 48.7 with respect to 46.1 °C of the constant interpolation case. Solving this problem takes about 2 minutes on a Pentium 4 3GHz, equipped with 2 GB of RAM. Increasing the number of variables could lead to better results because the interpolation would be more precise, however the optimization procedure could have some problems in reaching the convergence and the time for solving the problem could increase too much. II.4 Conclusion The problem of determining the temperature in a domain and along its boundary, using supplementary information along a part of the boundary of the domain, has been formulated as an optimization problem, minimizing a suitable functional. Convergence of the numerical solution has been obtained, starting from various guess solutions for the temperature along the boundary of the subdomain. The methodology can be usefully applied to predict the temperature of tumoral tissues during a thermal ablation treatment. In view of a real‐life application, the runtime of the algorithm should be optimized in terms of both hardware platform and software implementation. In fact, the treatment takes about ten‐twenty minutes and the reconstruction of the temperature should keep a computation time within a few seconds. 36
Chapter III CASE STUDY 2: SYNTHESIS OF COILS FOR MAGNETIC SIMULATION OF PERIPHERAL NERVES
37
III.1 Introduction Magnetic stimulation is a known method for stimulating excitable tissues with an electric current induced by an external time‐varying magnetic field (Malmivuo and Plonsey, 1993). It is important to note that both electric and magnetic stimulation methods excite the nervous membrane with electric current. The former does that directly, but the latter does it with the electric current which is induced within the conducting volume by the time‐varying applied magnetic field. The magnetic field penetrates unattenuated through regions such as the electrically insulating skull; this makes it possible to avoid a high density of stimulating current at the scalp in stimulating the central nervous system and thus avoid painful sensations; on the contrary the current injected by an electric stimulator, flowing directly through the skin, makes the stimulation painful. Moreover, no physical contact of the stimulating coil with the target tissue is required, unlike with electric stimulation. For these reasons, magnetic stimulation of peripheral nerves is used instead of the electrical stimulation for several clinical applications. Transcranial magnetic stimulation is currently under study as a treatment for severe depression, mania, auditory hallucinations (e.g., associated with schizophrenia), post‐traumatic stress disorder, obsessive‐compulsive disorder, generalized anxiety disorder and migraine headaches. It is also used to investigate the functionality of the different regions of the brain. However, the main use of this technique is as an action potential elicitor to investigate the conductivity of the nervous fiber. The piece of nervous fiber to be investigated is that which lays between the stimulation point (brain or spine/spinal roots or distal nerve) and the registration point (whatever muscle or nerve in the body). The commercially available devices for magnetic stimulation are characterized by a discharge capacitor connected to a resistor and to an inductive load (coil) in series with a thyristor. With the capacitor first charged at 2−3 kV, the gating of the thyristor into the conducting state will cause the discharging of the capacitor through the coil. The resulting current waveform can change from device to device, depending on the values of the capacitance, the resistance and the inductance of the circuit. In the so called “biphasic stimulator” the current 38
waveform is typically a damped sinusoidal pulse that lasts about 300 µs and has a peak value of 5−10 kA. The induced electric field sets free charges into coherent motion both in the intra‐ and extracellular spaces. Basically, any part of the cell membrane interrupting this motion of the charges becomes depolarised or hyperpolarised. In practice, however, the basic cellular mechanisms are unclear, although the macroscopic electromagnetic fields are well understood. The magnetic stimulation of the brain is due to the induced electric field in the nervous tissue (Malmivuo and Plonsey, 1993). The stimulation of peripheral nerves (in the case of infinite length and straight axons) is due to the so called “activating function” that is the derivative of the x component of the induced electric field with respect to x, if x is the direction in which the nervous fiber lays, at a prescribed time instant (Rattay, 1986; Roth and Basser, 1990). Some authors have found that the activating function as previously described is incomplete and they propose the addition of a term related to the field component transverse to the axon (Ruohonen, Panizza et al., 1996; Ruohonen, Ravazzani et al., 1996). However, in this work, the classical activating function from Rattay is taken into account (Rattay, 1986; Roth and Basser, 1990; Plonsey and Barr, 1995). Starting from some experiments, the idea of improving the existing device has evolved. Hence, the effects of the current distribution in the coil has been investigated. In the past, some authors investigated the coil shape to improve this kind of device (Hsu and Durand, 2001; Zimmermann and Simpson., 1996). However, the design procedure was based on repeated analyses, by varying a design variable at a time (Esselle and Stuchly, 1995; Hsu et al., 2003) while the magnetic analysis of the device relied on analytical models (Carbunaru and Durand, 2001; Davey and Epstein, 2000; Roth and Basser, 1990). Here, a procedure of automated design of the coil has been applied. Moreover, while in the past the attention was focused on the improvement of brain stimulation (Cohen et al., 1990; Kammer et al., 2001; Ravazzani et al., 2002), the solutions proposed in this work are suitable for the stimulation of peripheral nerves. 39
III.2 Experiments It is well known that coils, having different shapes, give rise to different effects on the tissue to be stimulated. In fact, different kinds of coils are commercially available (www.medtronic.com). As an example, the parabolic coil is used for the stimulation of the inner part of the brain, because it gives a deep stimulation; the butterfly coil is used to stimulate localized regions in the brain, because it allows a selective stimulation. Different shapes of coils give rise to different distributions in the space of the induced electric field and, consequentially, of the activating functions. In order to understand this behaviour a 3D finite element model (described in chapter III.3.2) was performed. Solving that model, it was possible to obtain the activating function for each coil modelled. The coils taken into account were: the circular, the parabolic and the butterfly (or figure of eight) coil, which are shown in Fig. III.1. Fig. III.1 – Circular, parabolic and butterfly coil. The isolevel curves of a typical activating function of a circular, a parabolic and a butterfly coil are shown in Figs. III.2, III.3 and III.4, respectively. These activating functions are calculated at the time instant when the current impulse starts, in a plane (x,z) parallel to the coil and 1 cm below it. 40
Fig. III.2 – Isolevel curves of a typical activating function of the circular coil. Fig. III.3 – Isolevel curves of a typical activating function of the parabolic coil. 41
Fig. III.4 – Isolevel curves of a typical activating function of the butterfly coil. Changing the coil, the site and the strength of the stimulation change. In order to test these effects, some experiments were made. III.2.1 The collision method revisited: procedure The ulnar nerve, a mixed motor and sensory nerve, is chosen for the stimulation. It was stimulated by means of the collision method, according to (Kimura, 1989). The traditional method consists of giving three electric impulses. In this method revisited, the second electric stimulation is replaced by a magnetic one (Alfonsi et al., 2005). The sites chosen for stimulations are the wrist and the axilla (a place between the elbow and the shoulder). In particular, a first stimulus is given at the wrist, the second and the third at the axilla. Because the second stimulus is magnetic, two different kinds of stimuli, a magnetic and an electric one, will be given at the axilla. Great attention is paid at positioning the coil in such a way that the sites of the two stimulations are the same. 42
The signals are recorded by surface electrodes. The active electrode is applied over the abductor digiti minimi muscle (ipotenar eminence) while the indifferent electrode is applied over the skin of the fifth digit (see Fig. III.5). Fig. III.5 – Experimental setup. When a stimulation is given, two action potentials are elicited: one will propagate in orthodromic direction, the other in antidromic direction. Thus, thanks to the first electric stimulation, the action potential that propagates in orthodromic direction (potential 1 in Fig. III.6) will be recordered by the surface 43
electrodes. The other action potential elicited (potential 1’) will collide with that elicited by the second stimulus applied at the axilla (potential 2’). These potentials will collide in a nerve segment located between the two stimulation points, as shown in Fig. III.6. Fig. III.6 – A first electric stimulation is given at the wrist, a second magnetic stimulation is given at the axilla. The coil is not represented. Successively, a third stimulus is applied at the axilla again, so that one of the potentials elicited (potential 3’) will pass through the nerve segment in which the collision took place (see Fig. III.7). Fig. III.7 – The collision takes place between the two stimulation points, the potential 3’ tests the fiber conditions. In this way, the third potential can test (from here the terms “test potential”) the fiber condition and in particular its refractoriness (Kandel and Schwartz, 2003). This experiment is repeated for each of the three coils available. 44
III.2.2 The collision method revisited: results The action potential is characterized by amplitude, duration and area. Fig. III.8 – Parameters characterizing an action potential. In order to compare the results, the parameter taken into account is the amplitude. However, the area is an integral parameter and hence stable as concerns possible measurement errors. It is also stable from a physiological point of view. The triplet of impulses is repeated different times: the time interval between the first and the second impulse is constant (4 ms in these experiments) to obtain complete collision in the nerve segment axilla‐wrist, while that between the second and the third stimulus is varied between 1 ms and 4.5 ms. When this time interval is too short, the third potential will not be able to pass through the region in which the collision occurred, because of the absolute refractoriness of the fibers; hence no signal due to the third stimulation will be registered at the 45
ipotenar eminence, as shown in Fig. III.9. Electrical artifact
Fig. III.9 – Signal recorded for ∆t = 1.5 ms: only the potential elicited by the wrist stimulation is visible. In Fig. III.9 an electrical artifact, that also occurs in Fig. III.10 and Fig. III.11, is highlighted. It is due to the stimulation electrode placed at the wrist, i.e. near the registration electrodes. Increasing the time interval between the second and the third stimulation, the test potential will pass through a fiber that from an absolute refractoriness period passes in a relative refractoriness period until the nerve fiber fully recovers its conductivity. Thus, a second action potential will be recorded and its amplitude will increase, increasing ∆t (see Fig. III.10 and Fig. III.11). Fig. III.10 – Signal recorded for ∆t = 2 ms: the potentials elicited by the wrist and by the axilla stimulation, respectively, are visible. 46
Fig. III.11 – Signal recorded for ∆t = 3.5 ms: the potential elicited by the axilla stimulation is visible and it has an amplitude greater than that for ∆t = 2 ms. The recovery (or excitability) curve is then reconstructed for each coil. It considers the amplitude of the second potential recorded with respect to the amplitude of a basal potential. The basal potential is a potential elicited by an electric stimulus, without any collision or other stimulations occurring before it. In Fig. III.12 the results are shown (Mognaschi et al., 2006). Fig. III.12 – Excitability curves of the three experiments made. Circular coil: blue square, butterfly coil: pink plus, parabolic coil: red circle. Considering one coil, e.g. the circular one, it can be noted that until ∆t=1.5 ms 47
no action potential, eventually elicited at the axilla, is recorded. This means that the fiber is in an absolute refractoriness period, because it does not allow the transit of the test action potential. From 1.5 to 2.5 ms, hence for 1 ms, the fiber is in a relative refractoriness period, because it allows the passage of an action potential, but the amplitude of it is reduced. From this experiment it seems that the circular and the butterfly coils have the same effect on the fiber, while the parabolic coil has the effect of prolungating both the absolute and the relative refractoriness periods. III.2.3 Measurements of the time constant of a nervous fiber with different transducers: procedure
For this experiment the stimulation point is set along the ulnar nerve, at the axilla, while the evoked potentials are recorded at the ipotenar of the same arm considered. Before starting the experiments, some parameters have to be measured. First, the minimal current amplitude, able to give the maximum amplitude of the evoked potential, is determined. A single electric stimulation is so given and the response of the nerve is recorded. The current is increased gradually and the nerve response is recorded at each time. When the amplitude of the curve recorded stops increasing, the maximal stimulation and the maximal evoked potential are determined. Second, some magnetic stimuli in order to measure the threshold of the nerve, are released to the nerve. Fixing the intensity of the stimulation, ten stimulations are repeated. The first stimulation has such an intensity that it is possible to record an action potential for each of the ten stimulations. The intensity of the stimulation is decreased gradually until it is possible to record less than six to ten action potentials. That value of the current amplitude is considered as the threshold of the nerve. The experiment consists of giving two stimulations: the first is magnetic (conditionating stimulus) and the second is electric (test stimulus). The first is under threshold (Plonsey and Barr, 2000) and has a time duration of 280 µs, while the second, delayed by ∆t with respect to the first, has an amplitude that 48
allows one to record an action potential of 50% of the maximal one, and it has a time duration of 70 µs. When ∆t tends to infinity, the second stimulus is not influenced by the first and the amplitude of the potential recorded is equal to 50% of the maximal one. Decreasing ∆t, the second stimulus is conditioned by the first (from here the term “conditioning stimulus”), in the sense that the fiber, first invaded by an electromagnetic field spread in a more or less wide region, changes its properties. As a consequence, the action potential recorded will have an increased amplitude: the conduction of the nerve results increased. For each delay ∆t considered, an action potential is recorded. Before and after this potential, a potential (basal potential) elicited by an electric stimulation without the conditioning stimulus is recorded. The conditioned action potential is compared with the average of the two basal potentials recorded before and after the conditioned potential. In order to avoid measurement errors and variability due to movements or physiological changes of the fiber, each action potential considered is the mean value of ten action potentials. This operation is directly made by the electromyograph device. In order to test the effect of the width of the region influenced by the electromagnetic field, the three coils mentioned in paragraph III.2 were used for this experiment. The coils were tested on two healthy male subjects, who were 24 (subject 1) and 18 (subject 2) years old. III.2.4 Measurements of the time constant of a nervous fiber with different transducers: results The parameter measured in the conditioned action potential (parameter C) is compared with the mean value of the same parameter measured in the previous (parameter Ap) and in the next ( parameter An) basal potential, for each ∆ti with i=1,…, n and n dependent on the experiment considered. In this way, a parameter per unit P is calculated: 49
Pi =
2C i
i = 1,..., n ( Ap i + An i )
(III.1) With reference to the amplitude as the parameter measured, the results obtained with the circular coil are shown in Fig. III.13: Fig. III.13 – Experimental (circle) and interpolated (continuous line) data for the circular coil. The experimental data are interpolated by an exponential law: p( ∆t ) = 1 + ke
− ( ∆t + T )
τ
(III.2) where T is the first time delay ∆t considered for the experiment and k and τ are identified by an automated procedure. In this way it is possible to calculate the time constant of the fiber τ. In literature the time constant was calculated using two electric stimuli (Bostock and Rothwell, 1997; Panizza et al., 1993). Moreover the procedure was a little different because the potential recorded always had the same amplitude. The 50
effect of the conductivity change of the nerve was evaluated by considering how much the current amplitude for the test stimulus had to be increased in order to record an action potential of the same amplitude. This time constant was calculated in the ways described previously, however, it is not the real time constant of the fiber; it is not an intrinsic parameter of the tissue. In fact it is called, more properly, strength duration time constant (Bostock and Rothwell, 1997). In this work, for simplicity, it will be called time constant. The results obtained for the two subjects examined are reported in Fig. III.14 e Fig. III.15: Fig. III.14 – Interpolated results referred to subject 1. Butterfly coil: blue line, parabolic coil: green line, circular coil: pink line. 51
Fig. III.15 ‐ Interpolated results referred to subject 2. Butterfly coil: blue line, circular coil: pink line. The parameters identified are: Table III.1‐ Parameters k and τ identified for the two subjects. Subject 1 Subject 2 Circular Parabolic Butterfly Circular Parabolic Butterfly coil coil coil coil coil coil k 0.13 0.16 0.17 0.38 ‐ 0.28 τ [µs]
208 102 35 40 ‐ 31 The order of magnitude of the time constant obtained by these experiments is the same as those measured by (Panizza et al., 1993). However a remark is necessary: both motor and sensory fibers belong to the ulnar nerve, but in these experiments the answer by the motor part of the nerve is evaluated. Because the time constant of the motor fibers is different from that of the sensory fibers (Panizza et al., 1993), a comparison has to be made with the time constant of the motor nerve. On the basis of this consideration, a time constant of 208 µs 52
(circular coil) seems too high for a motor fiber, while it seems to be a typical value for sensory fibers. The time constant measured for subject 2 with the circular coil seems to be more correct than that measured for subject 1. However, to draw any kind of conclusion, many more experiments and a greater number of subjects have to be considered. The conclusion that these two kinds of experiments suggests is an experimental confirmation that the shape of the coil has an important effect on the nervous fiber. Hence, from a medical point of view, the choice of the magnetic transducer can influence the result of clinical examinations and evaluations of pathologies. 53
III.3 Automatic procedure for the synthesis of the optimal coil The problem concerning magnetic stimulation is the space diffusion of the stimulus. This diffusion, as seen in the previous paragraphs, can lead to erroneous conclusions on the fiber properties. Moreover, it is responsible for a certain difficulty in stimulating a specific nerve, without stimulating another nerve that lies near the stimulated one. Hence, the characteristics that a good magnetic stimulator should have are a focalized and powerful stimulation. In this work, the performances of the circular coil are improved, acting on the current distribution of a cross‐section of the coil. Changing the shape of the coil will vary its inductance. This change will lead to a different current shape. The shape of the current, and in particular the derivative of the current with respect to time, is proportional to the induced electric field and, consequently, to the activating function. Changing the current flowing in the coil is not beyond the purposes of this work; hence, because the other components of the device (e.g. capacitors, resistors etc...) remain the same, the inductance of the coil must remain the same as that of the starting device (i.e. the circular coil). Thus, a constraint on the coil inductance is considered. In order to obtain robust solutions, a further constraint on the sensitivity of the objective functions is taken into account. III.3.1 The direct problem The direct problem to be solved is ∂A
− ∇ 2 A + µσ
(III.3) = µJ , in Ω , t > 0 ∂t
with A (x,y,z,t) magnetic vector potential, J (x,y,z,t) impressed current density,
µ magnetic permeability, σ electric conductivity and Ω a simply‐connected domain. Eq. (III.3) can be solved once the initial conditions and the boundary 54
conditions are assigned. The initial condition A (x,y,z,0) = 0 and the boundary conditions n × ( ∇ × A ) = 0 along ΓN (III.5) (III.6) n ⋅ ( ∇ × A ) = 0 along ΓD (III.4) are imposed. If the sources are localized in a finite region τ of the domain, which is considered to be infinite, and the domain is homogeneous, Eq. (III.3) has the solution µ
Jdτ
A=
(III.7) ∫
4 π τ r − rʹ
where r is the position of the point in which the vector potential is calculated,
rʹ the position of the volume element dτ. Knowing the vector potential A , the electric field E can be obtained by solving the following equation ∂A
E = −
(III.8) ∂t
from which it follows E=−
∂⎛ µ
Jdτ ⎞
⎜ ∫
⎟ ∂t ⎜⎝ 4 π τ r − rʹ ⎟⎠
(III.9) Hence, the three components of the electric field Ex(t), Ey(t) and Ez(t) can be 55
calculated Ex = −
∂ ⎛ µ J x dτ ⎞
⎜
⎟
∂t ⎜⎝ 4π ∫τ r − rʹ ⎟⎠
(III.10) Ey = −
∂ ⎛ µ J y dτ ⎞
⎜
⎟
∂t ⎜⎝ 4π ∫τ r − rʹ ⎟⎠
(III.11) Ez = −
∂ ⎛ µ J z dτ ⎞
⎜
⎟
∂t ⎜⎝ 4 π ∫τ r − rʹ ⎟⎠
(III.12) It can be noted that each component of the electric field exists if the corresponding component of the source J is different from zero. In the real device, the current flowing in the coil and the associated current density have the components x and z not null, while the y component is zero. Hence, the electric field (see Eqs. (III.10)‐(III.12)) has only the two components Ex and Ez different from zero in all the domain, except for x = 0, in which Ex = 0, and z = 0 in which Ez = 0, for symmetry properties of the problem considered. Let’s suppose a nerve laying in the x direction and a coil parallel to it, see Fig. III.16. Air
Coil
Forearm
Nerve
y
z
x
Fig. III.16 – 2D model of the coil and the arm. 56
In order to stimulate the nerve, it is necessary to have the activating function, and hence the electric field Ex, different from zero. The x component of the current density must be different from zero. Let’s suppose performing a 2D model of the coil and the nerve (see Fig. III.16). In particular, modelling the cross section of the coil in the x‐y plane with the current flowing in the z direction. The z‐component is the only component of the electric field not null (see Eq. (III.12)) and Ex and Ey are null in all the domain. In such a way, the stimulation of the nerve does not occur. On the contrary, modelling this problem with a 3D model, it is possible to assign two components (x and z) to the current source; so far, a nerve lying in the x direction results stimulated. Lets consider a nerve lying in the x direction with y= y and z= z and a 2D coil lying in the x‐z plane as shown in Fig. III.17. Coil
y
Coil z y φ
z
x
z
z x Nerve Nerve Fig. III.17 – Scheme of the position of the nerve and the coil. The Cartesian axes have the centre coincident with that of the coil and φ is the angle that the ray of the coil has with the direction (x,0,0), moving along the coil circumference. The induced electric field in x direction can be calculated as follows: 57
E x ( x) = −
µ 0 dJ 2 π
sin(φ)
dφ (III.13) ∫
2
4π dt 0 ( x − r cos( φ)) + ( y ) 2 + ( z − r sin(φ)) 2
This integral can be solved by integration per series and the solution exists. This problem is a typical three dimensional problem, hence a 2D model is not adequate. Once Ex is computed, the activating function Af can be calculated as: ∂E x
(III.14) Af =
∂x
Typical shapes for Ex and Af, calculated analitically in a plane x‐z, parallel to and 1 cm below a circular coil, are shown in Figs. III.18 and III.19. Fig. III.18 – Ex calculated analitically 1 cm below the coil plane at t=0+. The projection of the coil on the x‐z plane is represented by the circle; the current flows in a clockwise direction. 58
Fig. III.19 – Activating function calculated analitically 1 cm below the coil plane at t=0+. The projection of the coil on the x‐z plane is represented by the circle. A slice of Fig. III.19, for z= z , is considered. A typical shape for the activating function so obtained is represented in Fig. III.20. f1
80%f1
f2
Fig. III.20 – Slice obtained for z= z from Fig. III.19. 59
The parameters considered for the activating function, computed at t=0+, are the value at the peak (f1 in Fig. III.20) and the width (f2 in Fig. III.20), at the threshold of 80% of the peak value f1. From a physical viewpoint, f1 is related to the stimulation strength and f2 to its selectivity in space. III.3.2 The model A 3D finite element model of 64,000 nodes and 365,000 tetrahedrons was made. The arm is modelled by a cylinder made of three different materials: the skin, the fat and the muscle, as shown in Fig. III.21. Fig. III.21 – 3D finite element model. 60
The coil is modelled as stranded so that the effect of eddy currents in the coil is neglected. In order to identify the current carried by the coil, some measurements and simulations were performed. Because of technical problems in reaching the clamps of the coils with a probe for measuring the current directly, an indirect measurement was performed. The induced voltage was measured by the coil in a sensing coil. For this purpose a three turn sensing coil was built; the area of the coil was 1.58 cm2. This coil was connected to an oscilloscope that could record signals on a floppy disk. For each kind of coil, different measurements were made; in particular, for each coil, the power released was imposed to be 10%, 50%, 90% and 100 % of the maximum power released. Three coils were taken into account and, for each of them, different positions of the sensing coil with respect to the exciting coil were considered. In Fig. III.22 three coils considered with their sensing coil are represented. Fig. III.22 – Position of the sensing coil (white one) with respect to each kind of coil. Once the induced voltage V is measured (see Fig. III.23), from the Faraday‐
Neumann‐Lenz law we have d( ∫ B ⋅ dS)
dΦ
=− S
V=−
(III.15) dt
dt
61
and under the hypothesis of B constant all over the sensing coil surface, it is possible to obtain the magnetic induction field B (see Fig. III.24), integrating the induced voltage V in time and dividing it by the surface S of the sensing coil. Fig. III.23 – Induced voltage V measured (circular coil). 62
Fig. III.24 – Induction magnetic field obtained applying Eq. (III.12) to induced voltage V measured (circular coil). In Fig. III.23 and Fig. III.24 it can be noted that the period of the induced voltage and of the induction field is 280 µs. Because B is proportional to the current flowing in the coil, the current has a period of 280 µs. For identifying the peak value of the sinusoidal curve of the current, the experimental setup was simulated by the finite element method. For the circular coil and the parabolic coil an axial‐symmetric simulation was performed. For the butterfly coil a 3D simulation was performed. In Fig. III.25 the geometries used for simulations are shown. Fig. III.25 – Geometries of the three coils simulated; in each figure the sensing coil is visible. 63
In the FE model, the current in each coil is regulated in order to have, in the centre of the sensing coil surface, the maximum value of the magnetic field B obtained from the induced voltage measured for the same coil. The results obtained for the currents in each coil are shown in Table III.2. Table III.2 – Experimental and simulated data for each coil considered. The power of the magnetic stimulation was set to 100%. Coil No. of Induction field Measured Induction Estimated turns measured [T] period field current [A]
[µs] simulated [T] Circular 12 1.5 280 1.58 6000 Parabolic 14 0.85 310 0.87 6000 Butterfly 2×10 1.3 280 1.38 6000 The results obtained with the power set to 10%, 50% and 90% are coherent with results in Table III.2. Hence, the typical current carried by the circular coil used in clinical applications is i(t)= 6 10 3 sin(2.24 10 4 t) A for 0 < t < 280 µs , i(t)=0 elsewhere; the coil has 12 turns. This kind of coil was taken as a prototype to be improved. The magnetic and electrical properties of biological tissues varies with the frequency; for our purposes, a frequency of 3.57 kHz was taken into account. The values of electrical conductivity σ and permittivity ε, used for the biological tissues of the model, are reported in Table III.3 (Gabriel et al., 1996; niremf.ifac.cnr.it/tissprop/). Table III.3 – Electrical properties of biological tissues used in the model. Tissue Conductivity σ [Sm‐1] Relative permittivity εr Skin dry 0.0002 1135 Fat 0.0234 4562 Muscle 0.3344 78403 64
All tissues are assumed to be non‐magnetic (µ=µ0).
The Neumann boundary condition Eq. (III.5) is applied to the surfaces x=‐12.5 and x=12.5 cm, while the Dirichlet boundary condition Eq. (III.6) is applied to the surfaces y=‐10, y=10, z=‐10 and z=10 cm. The activating function is always evaluated for t=0+ and in the plane y=3.3 cm, i.e. 0.9 cm below the coil. After calculating for what value of z= z the peak exists, the functions f1 and f2 are evaluated on the slice obtained for z= z ; z varies with the coil shape. Because the activating function comes from a numerical derivation of the induced electric field, it is affected by a white noise, as shown in Fig. III.26. Fig. III.26 – A typical activating function affected by white noise. Fortunately, the frequency of the noise is much higher than that of the activating function; then, a butterworth filter of the fourth order is applied. After filtering, the activating function becomes that shown in Fig. III.27: 65
Fig. III.27 – The activating function of Fig. III.26 after the application of a digital filter. The parameters f1 and f2 (see Fig. III.20) will be expressed in Vm‐2 and cm. III.3.3 Inverse problem In order to synthetise the field source, a rectangular design region (see Fig. III.28) is considered and different distributions of current inside this region are taken into account. The design region has vertexes V1 = (1, 4.2), V2 = (5.5, 4.2), V3 = (5.5, 5.7) and V4 = (1, 5.7) cm. The current density and the linkage current are prescribed and equal to 266 Amm‐2 and 72000 At, respectively, as in the existing device, i.e. the circular coil. These constraints lead us to consider the winding area constant and, in particular, equal to 2.7 cm2. The upper contour of the distribution of current may exhibit either parabolic or linear shape. The distribution of current inside the design region is governed by four design variables [x1, x2, x3, x4] that are the coordinates of points P1 [x1, x2] and P2 [x3, x4] in Fig. III.28. In section, the sides of the coil are always vertical lines. 66
Design region
Air
V3
V4
P1
P2
Coil
V1
V2
Forearm
y
x
Fig. III.28 – Cross section of the coil. The design region is highlighted. The inverse problem can be formulated as follow: Formulation 1 ‐ acting on the design variables, find the distribution of current within the given design region such that the peak value of the activating function f1(x1, x2, x3, x4) (see Fig. III.20) is maximum and its width f2(x1, x2, x3, x4) (see Fig. III.20) is minimum at t=0+, under the constraints that the linkage current and the current density are constant. Successively, a second inverse problem is solved: Formulation 2 ‐ given the formulation 1, an additional constraint is added, such that the inductance of the coil must be equal to the inductance of the prototype within a ± 5% tolerance (Di Barba, Mognaschi and Savini, 2006). The inductance of a coil has to be calculated, normally, many times at each iteration. In order to limit the computational cost, the inductance of the coil is calculated by means of an axialsymmetric finite element model. In this way, the computational costs are comparable with a 2D finite element analysis. The prototype considered as the starting point has an inductance equal to 4.6 µH. Hence, the limits for the inductance constraint are [4.14, 5.05] µH. 67
A third inverse problem is then solved: Formulation 3 ‐ given the formulation 2, an additional constraint is added, such that the sensitivity calculated at each iteration must be lower than that calculated for the prototype. Traditionally, the evaluation of the sensitivity is very expensive, in terms of objective function calls (Haslinger and Neittaanmäki, 1988). Here a free‐extra cost for the evaluation of the sensitivity is described (Di Barba and Mognaschi, 2006). Given a set of np>>1 points sampling the nv‐dimensional design space, a cluster is defined as the subset formed by the points such that their distances d i,j
⎡ ⎛
⎢ n v ⎜ x i ,k − x j ,k
= ⎢∑⎜
k =1 sup x − inf x
k
⎢ ⎜⎝ i k
i
⎣
⎞
⎟
⎟
⎟
⎠
2
⎤
⎥
⎥
⎥
⎦
1
2
, i = 1, n p − 1 , j = 2 , n p , j > i (III.16) are lower than the mean radius ρ defined as n p −1 n p
2d i , j
ρ= ∑ ∑
i =1 j= i + 1 n p ( n p − 1)
(III.17) Given the current four‐dimensional design vector ξ , the cluster, which ξ belongs to, is considered by means of (III.16) and (III.17). The sensitivity is then evaluated according to the ∞‐norm: (III.18) s λ (ξ ) = sup sup fi (ξ(λ, j)) − fi (ξ(λ, k )) i =1 , n f
j =1 , n s
k =1 , n s
j≠ k
where ℓ is the iteration index, nf=2 is the number of objectives and ns>1 is the number of cluster elements. 68
In this way, at each iteration the sensitivities of both the objectives are evaluated on the basis of the history of the optimization. The biggest value, normalised with respect to the starting point, is taken as the sensitivity of the current design. The approach described is inexpensive in the frame of an optimization process; indeed, the bigger the number of cluster elements, the more accurate the sensitivity. Nonetheless, the evaluation might be not applicable during the first iterations, when the cluster includes just the current element. In particular, the sensitivity of the starting point is calculated beforehand after perturbing the initial design vector. To this end, a cluster of e.g. ns=10 individuals is randomly generated in the feasible region and, on the basis of this cluster, the sensitivity is calculated according to (III.18). III.3.4 Numerical results The starting point of all the minimizations is the circular coil, the cross section of which is defined by x0=[1, 4.8, 5.5, 4.8] cm. The activating function is represented in Fig. III.29: 69
Fig. III.29 – Activating function calculated for the circular coil. The objective functions for the starting point are f0=[3749, 2.76]. Problem 1: The arrival point for the first problem is x=[2.47, 5.54, 5.19, 4.9] and f=[4817, 2.15]. The objective functions f1 and f2 are improved by 28.5% and 22.1%, respectively. The shape of the coil is shown in Fig. III.30: 70
Fig. III.30 – Shape of the optimal coil for problem 1. The activating function is shown in Fig. III.31. For symmetry reasons, the only positive half‐wave is represented. Starting geometry
Arrival geometry
Fig. III.31 – Activating function of the starting geometry and of the arrival geometry for problem 1. The cross section of each coil is represented. 71
The history of the objective functions is shown in Fig. III.32 in which the functions are normalized with respect to the value of the starting point: Fig. III.32 – History of the objective functions, f1: blue line, f2: pink line. Problem 2: The arrival point for the second problem is x=[1.44, 5.39, 4.51, 5.53] and f=[3907, 2.1]. The objective functions f1 and f2 are improved by 4.2% and 23.9%, respectively. The inductance of this coil is 4.86 µH. The shape of the coil is shown in Fig. III.33: 72
Fig. III.33 – Shape of the optimal coil for problem 2. The activating function is shown in Fig. III.34: Starting geometry
Arrival geometry
Fig. III.34 – Activating function of the starting geometry and of the arrival geometry for problem 2. The cross section of each coil is represented. 73
The history of the inductance is represented in Fig. III.35. Fig. III.35 – History of the inductance. For all the iterations, the inductance is within the limits, i.e. [4.14, 5.05] µH. The history of the objective functions is shown in Fig. III.36 in which the functions are normalized with respect to the value of the starting point: 74
Fig. III.36 – History of the objective functions, f1: blue line, f2: pink line. Until the fiftieth iteration, the dispersion of the method does not vary because the optimization method stores an optimization history. In these iterations it tests the possible directions. In Fig. III.36 it can be noted that the trend, from ten to fifty iterations, is flat. This is due to the inductance constraint. In fact, for this problem, it is a tight constraint and the algorithm is able to choose only the same point for a lot of iterations. In fact, for problem 1 (see Fig. III.32) this behaviour cannot be noted. Problem 3: For this problem the inductance constraint is observed at each iteration, as can be seen in Fig. III.37. 75
Fig. III.37 – History of the inductance. The sensitivity calculated at each iteration is shown in Fig. III.38. Fig. III.38 – Sensitivity calculated at each iteration. It can be noted that at each iteration the maximum value of the sensitivity is 76
under the value of the initial sensitivity, i.e. 0.13. The arrival point for the third problem is x=[1.82, 5.57, 4.31, 4.54] and f=[3969, 2.13]. The objective functions f1 and f2 are improved by 5.9% and 22.8%, respectively. The inductance of this coil is 4.62 µH.
The shape of the coil is shown in Fig. III.39: Fig. III.39 – Shape of the optimal coil for problem 3. The activating function is shown in Fig. III.40: 77
Starting geometry
Arrival geometry
Fig. III.40 – Activating function of the starting geometry and of the arrival geometry for problem 3. The cross section of each coil is represented. The history of the objective functions is shown in Fig. III.41 in which the functions are normalized with respect to the value of the starting point: Fig. III.41 – History of the objective functions, f1: blue line, f2: pink line. 78
The objective function space was investigated. A random sampling was performed. Fig. III.42 – Sampling of the objective space (dot), random sampling with inductance constraint (circle), starting point (pink square), arrival point for problem 1 (green asterisk), arrival point for problem 2 (red triangle) , arrival point for problem 3 (blue cross). It can be noted that in Fig. III.42 the objective functions f1 and f2 can vary by about 50% and 30%, respectively. The points with the inductance constraint seem to stay in a well‐defined subregion of the objective space. This subregion 79
has a variation for f1 comparable with that of the objective space sampled, while for f2 the possible variation is limited (about 20%). The arrival point for problem 1 improves both the objective functions. The solution of problem 2, because of the inductance constraint, does not improve function f1 a lot, while it improves function f2. The solution of problem 3 is close to that of problem 2, but this solution is robust, because it fulfils the sensitivity constraint. III.4 Conclusion A number of experiments regarding the magnetic stimulation of the ulnar nerve in vivo were made. These experiments confirmed what is shown in literature i.e. the shape of the coil influences the stimulation. From this idea, an automated procedure to improve the design of a magnetic stimulator was developed and applied to a commercially available device. The coil producing the optimal activating function was determined by the implemented optimization procedure. Then, a further inductance constraint was added to the optimization procedure and the optimal coil for this problem was found. Finally, a robust optimal coil was obtained by adding a sensitivity constraint to the optimization problem; the calculation of the sensitivity is free of extra costs. The objective space was also explored in order to compare the different solutions. 80
Conclusion 81
Two different kinds of inverse problems in bioelectromagnetism have been studied. The first belongs to the identification problems. The problem was to reconstruct the temperature on a slice of a patient torso, during a thermoablation treatment of liver tumor. The known data are the temperatures measured by a certain number of thermocouples placed on the patient’s skin. From a CT image, the domain of the problem was built. The direct problem was solved by the finite element method. First, the effect of capillaries, uniformly distributed all over the domain, was considered to take into account the perfusion phenomenon. The boundary of the tumor was considered to be isothermal; despite the roughness of this assumption, good results were obtained. This first case was useful in testing the identification algorithm. A second case, in which the localized perfusion was also considered, was solved: two big vessels close to the tumor were considered in the domain. This more realistic case did not give us good results as the simpler case, considering the boundary of the tumor to be isothermal. Hence, a bilinear approximation of the variable temperature along the tumor boundary was considered. In this way, the error in reconstructing the temperature was acceptable. Increasing the variables of optimization, the accuracy would increase. However, the optimizer would not be able to work with such a large number of variables. Using three variables, the computational costs are quite compatible with a thermal ablation treatement. In order to apply this method to a real thermal ablation treatment, first of all it is necessary to improve the computational cost: the time consumed for the prediction must be much shorter than the time used for the ablation. Then, the accuracy must be improved. However, because the range of temperature useful for the thermal ablation therapy is large and it is important for the temperature to be over 42 °C, sub‐estimating the temperature, allows us to reach correct temperatures inside the tumor, without reaching 100‐105 °C in which there would be bubble formation. The second case study considered came from some experiments made on the peripheral nervous system. The magnetic stimulation of the ulnar nerve was performed, using different stimulation methods. The first was the collision method and the second a method to measure the time constant of the fiber. For 82
each method, three different coils were tested. 3D finite element models of the coils showed the different effect of each transducer on the nervous fiber. This effect was due to the shape of the coil considered: changing the coil shape, the stimulation changes. In order to have a strong and selective stimulation, a design of the optimal coil, acting on the coil shape, was performed. The automated procedure was based on a 3D finite element model. Three coil designs were performed: the first used an unconstrained optimization to minimize the two objective functions. In the second, a constraint on the inductance of the coil was added and in the third a further constraint on the sensitivity was added. From these optimizations, three different coil shapes emerged, each of them performing more or less the objective functions. In particular the coil designed with an unconstrained optimization improved the objective functions by 29 and 22 %, respectively. The coil with the inductance constraint substantially improved only the selectiveness of the stimulus and not the strength very much. Also the coil on which a sensitivity analysis was performed, improved the selectiveness of the stimulus and not the strength too much. However, the latter implements a robust solution, due to the sensitivity constraint. In this work, methodologies which have been traditionally applied to classical problems of electromagnetism, were presented and developed for solving two classes of inverse problems originating in bioelectromagnetism. The results obtained were satisfactory and encouraging; in particular, original, unpredicted solutions have been found. 83
Bibliography 84
Alfonsi, E., Di Barba, P., Mognaschi, M.E., and Savini, A., (2005), Effects of pulsed magnetic fields on motor fibers of a peripheral nerve: experimental results and simulations, Proceedings of XII International Symposium on Interdisciplinary Electromagnetic, Mechanic & Biomedical Problems ISEM 2005, Bad Gastein (Austria). Bostock, H., and Rothwell, J.C., (1997), Latent addition in motor and sensory fibres of human peripheral nerve, Journal of Physiology, Vol. 488.1. Brauer, H., Haueisen, J., Ziółkowski, M., Tenner, U., and Nowak, H., (2000), Reconstruction of extended current sources in a human body phantom applying biomagnetic measuring techniques, IEEE Transactions on Magnetics, Vol. 36, No. 4. Carbunaru, R., and Durand, D.M., (2001), Toroidal coil models for transcutaneous magnetic stimulation of nerves, IEEE Transactions on Biomedical Engineering, Vol. 48, No. 4. Carcangiu, S., Di Barba, P., Fanni, A., Mognaschi, M.E., and Montisci, A., (2006), Comparison of multi‐objective optimisation approaches for inverse magnestostatic problems, in press on COMPEL (The International Journal for Computation and Mathematics in Electrical and Electronic Engineering). Cavaliere, V., Formisano, A., Martone, R., and Primizia, M., (2000), A genetic algorithm approach to the design of split coil magnets for MRI, IEEE Transactions on Applied Superconductivity, Vol. 10, No. 1. Cohen L.G., Roth, B.J., Nilsson, J., Dang, N., Panizza, M., Bandinelli, S., Friauf, W., and Hallett, M., (1990), Effects of coil design on delivery of focal magnetic stimulation. Technical considerations, Electroencephalography and Clinical Neurophysiology, Vol. 75. 85
Colli Franzone, P., (1980), Regularization Method Applied to an Inverse Problem in Electrocardiology, Computing Methods in Applied Sciences and Engineering. D’Ambrosio, V., Di Barba, P., Dughiero, F., Mognaschi, M.E., and Savini, A., (2005), Non‐invasive thermometry for the thermal ablation of liver tumor: a computational methodology, Proceedings of XII International Symposium on Interdisciplinary Electromagnetic, Mechanic & Biomedical Problems ISEM 2005, Bad Gastein (Austria) . D’Ambrosio, V., Di Barba, P., Dughiero, F., Mognaschi, M.E., and Savini, A., (2006), Identification of boundary conditions in a problem of thermal ablation of liver tumor, Proceedings of IX International Workshop on Optimisation and Inverse Problems in Electromagnetism OIPE 2006, Sorrento. Davey, K., and Epstein, C.M., (2000), Magnetic stimulation coil and circuit design, IEEE Transactions on Biomedical Engineering, Vol. 47, No.11. Di Barba, P., (2001), A fast evolutionary method for identifying non‐inferior solutions in multicriteria shape optimisation of a shielded reactor, COMPEL (The International Journal for Computation and Mathematics in Electrical and Electronic Engineering), vol. 20, no. 3. Di Barba, P., and Mognaschi, M.E., (2004), Recent experiences of multiobjective optimisation in electromagnetics: a comparison of methods, COMPEL (The International Journal for Computation and Mathematics in Electrical and Electronic Engineering), Vol. 24, No. 3. Di Barba, P., and Mognaschi, M.E., (2006), Optimal shape design of a coil for nerve stimulation: a sensitivity based approach, Proceedings of XII International IGTE Symposium on Numerical Field Calculation in Electrical Engineering IGTE 2006. 86
Di Barba, P., Mognaschi, M.E., and Savini, A., (2006), Synthesis of a field source for the electromagnetic stimulation of peripheral nerves, Proceedings of IX International Workshop on Optimisation and Inverse Problems in Electromagnetism OIPE 2006, Sorrento. Dughiero, F., and D’Ambrosio, V., (2004), FEM Models of Radiofrequency Thermal Treatments in Cancer Cure, Proceedings of Conf. HES 04 – Heating by Electromagnetic Sources, Padova. Esselle, K.P., and Stuchly, M.A., (1995), Cylindrical tissue model for magnetic field stimulation of neurons: effects of coil geometry, IEEE Transactions on Biomedical Engineering, Vol. 42, No. 9. Gabriel, S., Lau, R.W., and Gabriel, C., (1996), The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz, Physics in Medicine and Biology, Vol. 41. Hadamard, J., (1923), Lectures on the cauchy’s problem in linear partial differential equations, Yale University Press. Haslinger, J., and Neittaanmäki, P., (1988), Finite element approximation for optimal shape design: theory and applications, Wiley. Hokland, S.L., Pedersen, M., Salomir, R., Quesson, B., Stødkilde‐Jørgensen, H., and Moonen, C.T.W., (2006), MRI‐Guided Focused Ultrasound: Methodology and Applications, IEEE Transactions on Medical Imaging, Vol. 25, No. 6. Hsu, K.H., and Durand, D.M., (2001), A 3‐D differential coil design for localized magnetic stimulation, IEEE Transactions on Biomedical Engineering, Vol. 48, No. 10. 87
Hsu, K.H., Nagarajan, S.S., and Durand, D.M., (2003), Analysis of Efficiency of Magnetic Stimulation, IEEE Transactions on Biomedical Engineering, Vol. 50, No. 11. Kammer, T., Beck, S., Thielscher, A., Laubis‐Herrmann, U., and Topka, H., (2001), Motor threshold in humans: a transcranial magnetic stimulation study comparing different pulse waveforms, current directions and stimulator types, Clinical Neurophysiology, Vol. 112. Kandel E. R., Schwartz J. H., (2003), Principi di Neuroscienze ‐ Terza edizione, Casa Editrice Ambrosiana, Gruppo Editoriale Zanichelli, Milano [English version: Principal of Neural Science, 1985, Elsevier Science Publishing CO., New York]. Kimura, J., (1989), Electrodiagnosis in Diseases of Nerve and Muscle: Principles and Practice, Edition 2, F.A. Davis Company – Philadelphia. Küçükaltun‐Yıldırım, E., Pantazis, D., and Leahy, R.M., (2006), Task‐based comparison of inverse methods in magnetoencephalography, IEEE Transactions on Biomedical Engineering, Vol. 53, No. 9. Lang, J., Erdmann, B., Seebass, M., (1999), Impact of nonlinear heat transfer on temperature control in regional hyperthermia, IEEE Transactions on Biomedical Engineering, Vol. 46, No. 9. Malmivuo, J., and Plonsey, R., (1993), Bioelectromagnetism, Butler. Online: http://butler.cc.tut.fi/~malmivuo/bem/bembook/index.htm Mognaschi, M.E., Alfonsi, E., Nilsson, J., and Moglia, A., (2006), Experimental results on excitability of motor axons of a peripheral nerve by means of a magnetic stimulus, Functional Neurology, Vol. XXI, No 2. 88
Neittaanmäki, P., Rudnicki, M., and Savini, A., (1996), Inverse problem and optimal design in electricity and magnetism, Oxford Science Publication ‐ Oxford University Press. Ngale Haulin, E., and Vinet, R., (2003), Multiobjective optimization of hand prosthesis mechanisms, Mechanism and Machine Theory, Vol 38. Pałka, R., (1983), Application of the Finite Element Technique to Continuation Problems of Stationary Fields, IEEE Transactions on Magnetics, Vol. Mag‐19, No. 6. Panizza, M., Nilsson, J., Roth, B.J., Rothwell, J., and Hallet, M., (1993), The time constants of motor and sensory peripheral nerve fibers measured with the method of latent addition, Electroencephalography and Clinical Neurophysiology, Vol. 93. Paruch, M., and Majchrzak, E., (2004), Identification of the Tumor Region on the Basis of Skin Surface Temperature, Proceedings of ECCOMAS 2004 – European Congress on Computational Methods in Applied Sciences and Engineering, Jyväskylä. Pennes, H., (1948), Analysis of tissue and arterial blood temperatures in the resting human forearm, Journal of Applied Physiology, Vol. I, No. 2. Plonsey, R., and Barr, C., (1995), Electric field stimulation of excitable tissue, IEEE Transactions on Biomedical Engineering, Vol. 42, No. 4. Plonsey, R., and Barr, C., (2000), Bioelectricity, a quantitative approach ‐ Second edition, Kluwer academic/Plenum Publishers. Preston, T.W., Reece, A.B.J., (1982), Solution of 3‐dimensional eddy current problems: the T‐Ω method, IEEE Transactions on Magnetics, Vol. Mag‐18, No. 2. 89
Rattay, F., (1986), Analysis of models for external stimulation of axons. IEEE Transactions on Biomedical Engineering, Vol. 33. Ravazzani, P., Ruohonen, J., Tognola, G., Anfosso, F., Ollikainen, M., Ilmoniemi, R.J., and Grandori, F., (2002), Frequency‐related effects in the optimization of coils for the magnetic stimulation of the nervous system, IEEE Transaction on Biomedical Engineering, Vol. 49, No. 5. Rieke, V., Ross, A.B., Nau, W.H., Diederich, C.J., Sommer, G., and Butts, K., (2004), MRI‐Temperature Mapping During Ultrasound Prostate Ablation Using Fat for Phase Estimation, Proceedings of the 26th Annual International Conference of the IEEE EMBS, San Francisco, CA, USA. Roth, B.J., and Basser, P.J., (1990), A Model of the Stimulation of a Nerve Fiber by Electromagnetic Induction, IEEE Transactions on Biomedical Engineering, Vol. 37, No. 6. Ruohonen, J., Panizza, M., Nilsson, J., Ravazzani, P., Grandori, F., and Tognola, G., (1996), Transverse‐field activation mechanism in magnetic stimulation of peripheral nerves, Electroencephalography and Clinical Neurophysiology, Vol. 101. Ruohonen, J., Ravazzani, P., Nilsson, J., Panizza, M., Grandori, F., and Tignola, G., (1996), A volume‐conduction analysis of magnetic stimulation of peripheral nerves. IEEE Transactions on Biomedical Engineering 43, 669–678. Tungjitkusolmun, S., Staelin, S.T., Haemmerich, D., Tsai, J.Z., Cao, H., Webster, J.G., Lee, F.T., Mahavi, D.M., and Vorperian, V.R., (2002), Three‐dimensional finite‐element analysis for radio‐frequency hepatic tumor ablation, IEEE Transactions on Biomedical Engineering, Vol. 49, No. 1. 90
Ulucakli, M.E., (2006), Simulation of Radiofrequency Ablation and Thermal Damage to Tissue, Proceedings of the 32nd Annual Northeast Bioengineering Conference. Zimmermann, K.P., and Simpson, R.K., (1996), Slinky coils for neuromagnetic stimulation, Electroencephalography and Clinical Neurophysiology, Vol. 101. Ziółkowski, M., and Brauer, H., (1996), Methods of mesh generation for biomagnetic problems, IEEE Transactions on Magnetics, Vol. 32, No. 3.
91
Papers authored by the candidate 92
Alfonsi, E., Di Barba, P., Mognaschi, M.E., and Savini, A., (2005), Effects of pulsed magnetic fields on motor fibers of a peripheral nerve: experimental results and simulations, Proceedings of XII International Symposium on Interdisciplinary Electromagnetic, Mechanic & Biomedical Problems ISEM 2005, Bad Gastein (Austria). Carcangiu, S., Di Barba, P., Fanni, A., Mognaschi, M.E., and Montisci, A., (2006), Comparison of multi‐objective optimisation approaches for inverse magnestostatic problems, Proceedings of IX International Workshop on Optimisation and Inverse Problems in Electromagnetism OIPE 2006, Sorrento. In press on COMPEL (The International Journal for Computation and Mathematics in Electrical and Electronic Engineering). D’Ambrosio, V., Di Barba, P., Dughiero, F., Mognaschi, M.E., and Savini, A., (2005), Non‐invasive thermometry for the thermal ablation of liver tumor: a computational methodology, Proceedings of XII International Symposium on Interdisciplinary Electromagnetic, Mechanic & Biomedical Problems ISEM 2005, Bad Gastein (Austria) . In press on IJAEM (International Journal of Applied Electromagnetics and Mechanics) D’Ambrosio, V., Di Barba, P., Dughiero, F., Mognaschi, M.E., and Savini, A., (2006), Identification of boundary conditions in a problem of thermal ablation of liver tumor, Proceedings of IX International Workshop on Optimisation and Inverse Problems in Electromagnetism OIPE 2006, Sorrento. Di Barba, P., and Mognaschi, M.E., (2004), Recent experiences of multiobjective optimisation in electromagnetics: a comparison of methods, COMPEL(The international journal for computation and mathematics in electrical and electronic engineering), Vol. 24, No. 3. Di Barba, P., and Mognaschi, M.E., (2006), Optimal shape design of a coil for nerve stimulation: a sensitivity based approach, Proceedings of XII International IGTE Symposium on Numerical Field Calculation in Electrical 93
Engineering IGTE 2006. Di Barba, P., Mognaschi, E.R., Mognaschi, M.E., and Savini, A., (2004), Identifying material properties of a dielectric motor, COMPEL (The international journal for computation and mathematics in electrical and electronic engineering) Vol.24, No. 3. Di Barba, P., Mognaschi, E.R., Mognaschi, M.E., and Savini, A., (2005), On a class of small dielectric motors: experimental results and shape design, A. Krawczyk, S. Wiak, and X.M. Lopez‐Fernandez, ʺElectromagnetic Fields in Mechatronics, Electrical and Electronic Engineeringʺ, IOS Press. Di Barba, P., Mognaschi, M.E., and Savini, A., (2004), Effetto di campi magnetici impulsati sulla refrattarietà della fibra nervosa umana, XX Riunione Nazionale dei Ricercatori di Elettrotecnica – ET2004, Salerno. Di Barba, P., Mognaschi, M.E., and Savini, A., (2005), Studio sperimentale sulla refrattarietà della fibra nervosa indotta da campi magnetici impulsati, XXI Riunione Nazionale dei Ricercatori di Elettrotecnica – ET2005, Roma. Di Barba, P., Mognaschi, M.E., and Savini, A., (2006), Synthesis of a field source for the electromagnetic stimulation of peripheral nerves, Proceedings of IX International Workshop on Optimisation and Inverse Problems in Electromagnetism OIPE 2006, Sorrento. Di Barba, P., Mognaschi, M.E., Savini, A., and Alfonsi, E., (2006), Evidenze sperimentali dell’effetto di forma del trasduttore sulla neurostimolazione magnetica, XXII Riunione Nazionale dei Ricercatori di Elettrotecnica – ET2006, Torino. Mognaschi, E.R., and Mognaschi, M.E., (2003), Earthquake source localisation by means of electromagnetic precursors, I International Workshop on Earthquake Prediction, Athens. 94
Mognaschi, M.E., Alfonsi, E., Nilsson, J., and Moglia, A., (2006), Experimental results on excitability of motor axons of a peripheral nerve by means of a magnetic stimulus, Functional Neurology, Vol. XXI, No 2. 95