Magnetic field effects on the vestibular system: calculation of the pressure on the cupula due to ionic current-induced Lorentz force

Large static magnetic fields may be employed in magnetic resonance imaging (MRI). At high magnetic field strengths (usually from about 3 T and above) it is possible for humans to perceive a number of effects. One such effect is mild vertigo. Recently, Roberts et al (2011 Current Biology 21 1635–40) proposed a Lorentz-force mechanism resulting from the ionic currents occurring naturally in the endolymph of the vestibular system. In the present work a more detailed calculation of the forces and resulting pressures in the vestibular system is carried out using a numerical model. Firstly, realistic 3D finite element conductivity and fluid maps of the utricle and a single semi-circular canal containing the current sources (dark cells) and sinks (hair cells) of the utricle and ampulla were constructed. Secondly, the electrical current densities in the fluid are calculated. Thirdly, the developed Lorentz force is used directly in the Navier–Stokes equation and the trans-cupular pressure is computed. Since the driving force field is relatively large in comparison with the advective acceleration, we demonstrate that it is possible to perform an approximation in the Navier–Stokes equations that reduces the problem to solving a simpler Poisson equation. This simplification allows rapid and easy calculation for many different directions of applied magnetic field. At 7 T a maximum cupula pressure difference of 1.6 mPa was calculated for the combined ampullar (0.7 µA) and utricular (3.31 µA) distributed current sources, assuming a hair-cell resting current of 100 pA per unit. These pressure values are up to an order of magnitude lower than those proposed by Roberts et al using a simplistic model and calculation, and are in good agreement with the estimated pressure values for nystagmus velocities in caloric experiments. This modeling work supports the hypothesis that the Lorentz force mechanism is a significant contributor to the perception of magnetic field induced vertigo.


Introduction
Magnetic field induced vertigo (MFIV) has been consistently reported by a number of staff working near and in strong magnetic fields (Schenck 1992) as used in magnetic resonance imaging (MRI). The vestibular system (VS) is known to play a crucial part in MFIV, evidence being that individuals (or animals) lacking labyrinthine function report none of the characteristic symptoms (Houpt and Houpt 2010). A number of possible mechanisms to explain VS stimulation have been enumerated before and these include: magnetic susceptibility of the vestibular maculae; induced currents causing galvanic vestibular stimulation (IGVS); as well as magneto-hydrodynamic (MHD) mechanisms (Glover et al 2007). However, experiments to determine the exact mechanism are not straightforward and suffer from subjective bias and the difficulty of performing good control and sham experiments. It is important to identify the precise principal mechanism (although it may be acknowledged that more than one may operate independently and simultaneously) as work-place legislation and regulations are formulated on the basis of likely acute effects on the human (ICNIRP (1998), HPA (2008), EU (2004) directives). It is certainly true from the above mechanisms that the VS is not working outside its limits or is damaged by the application of the magnetic field. Any vertigo or dizziness sensations are likely to be the result of vestibular inputs to the brain rather than a direct consequence of magnetic fields interacting with the central nervous system.
The recent findings by Roberts et al (2011) have introduced a new perspective to this problem. They consider the Lorentzian force generated by the interaction of the static magnetic field and the ionic (potassium flux) currents present in the endolymph which subsequently produce a pressure on the cupula. They are able to calculate a figure for the pressure exerted on this membrane relating not to the magnetic gradient field or the temporally changing field, but to the static magnetic field. Some simple approximations are made in their calculation, the most important being that the current flows purely linearly across the ampulla and hence the Lorentz force has a purely orthogonal orientation to the cupula. However, this geometry and assumption is not realistic and the current paths due to potassium cycling through the dark cells (current sources) and hair cells (current sinks) are much more complex. It is also possible to argue that much of the Lorentz force simply cancels without producing a net pressure. Hence it is imperative that a more realistic model now be constructed in order to calculate the current densities, Lorentz forces in the endolymph and resultant cupula pressures. This work achieves this aim by designing a realistic model and performing an exact numerical calculation, providing for the first time a more realistic estimate of the values of the pressure differential across the cupula due to the Lorentz mechanism.

Theory
Two independent but geometrically related models were used to run electromagnetic and fluid dynamics simulations. The calculation is carried out by employing three sequential steps.
(1) The calculation of the current density field using a quasi-static electromagnetic simulation.
This uses a conductivity map and places current sources and sinks in order to calculate the scalar potential field, and from it, the current density field. (2) The Lorentz force field calculation based on a specific orientation of a static uniform magnetic field. (3) The fluid dynamics calculation based on the Navier-Stokes equation, using the geometry of the problem and the Lorentz force field.

Quasi-static simulation
The Maxwell equations for a static magnetic field (B), in the quasi-static approximation hold: where E is the electric field, J is the current density, ρ is the charge density and ε 0 is the permittivity of free space. In a conductive fluid such as the endolymph there are no free charges and any charge accumulation at cell walls is ignored (ρ = 0). In this model the current sources are too small to generate a magnetic field which perturbs the applied field and hence ∇ × B = 0. For equation (3), there exists a scalar potential φ such that, Under steady state conditions, and knowing that there is no net charge density, ∇ · E = 0, leading to a solenoidal current density ∇ · J = 0. The distributed source of static current density can be defined by punctual sources that are added to the equation, such that where s is the volume distribution of current sources, in units of A m −3 . Since the current density is defined as where σ refers to the electrical conductivity with units of S m −1 . Taking the combinations of equations (5)-(7) and applying the divergence theorem leads to the integral form, where S is a closed surface, dA refers to the differential surface element, and V the volume enclosed by the surface S. Discretization of this equation in isotropic, iso-volumetric voxels leads to a linear system of equations in the standard form Mx = b, where M is a sparse matrix containing the conductivity and connectivity information for each of the element's surfaces, b contains the source current distribution and magnitudes, and x is the potential φ. A finite volume method approach for the volume discretization was adopted, ensuring the continuity equation is preserved for each voxel, and that each element only connects to its six nearest neighbors, turning M into a very sparse matrix, with all elements stored in seven diagonals. Once solved, current density can be found by taking J = −σ ∇φ. The BiCGSTAB algorithm (Van der Vorst 1992) was used to solve the linear system, using as the stopping criterion for convergence, a residual less than 10 −10 . The algorithm execution is highly parallelizable, so a C++/CUDA code was written to take advantage of this feature. CUDA is a computing architecture that can use Nvidia graphical processor units (GPU) cards to perform calculations in parallel (Nvidia Corporation, Santa Clara, CA USA). The simulations were tested and compared against different methods and analytical solutions: potentials on a square plate (van der Pauw 1958); the electric field on the surface of a sphere due to a transverse field (Bencsik et al 2002); and the potential on the surface of nested spheres due to a dipole (Zhang 1995).

Fluid dynamics equations
The Navier-Stokes equations for fluids provides the most general description of the dynamics of compressible, viscid fluids, where v is the fluid velocity, ρ is the fluid density (kg m −3 ), p is the pressure, f is the force density (N m −3 ) and μ is the fluid viscosity. Endolymph is essentially an incompressible fluid with low viscosity (μ 1) and we are interested in finding the steady state solution, so the equations can be reduced to ∇ p = f − ρv · ∇v + μ∇ 2 v.
(10) As an additional approximation we propose to consider the force per unit volume f to be considerably larger than the advective acceleration v · ∇v. Since f is proportional to J (∼1 A m −2 ) and B (7 T), and maximum advective acceleration is ∼10 −10 m s −2 , the approximation seems reasonable. The approximation breaks for regions where J and B are almost parallel, but in this situation pressures are low and the relative error introduced is of less consequence. This leads to (11), with solution similar to the electromagnetic case.
(11) Simulations with the complete equations were also carried out and compared to the approximated equations in order to justify this approach, with results for the cupular pressure diverging no more than 15% between methods. This level of divergence is deemed acceptable, particularly as the model geometry has a number of assumptions and even unknowns as discussed in the next section. The discretization scheme is now analogous to the electromagnetic equations, leading to a similar linear system of equations, which can also be solved for the pressures with the BiCGSTAB algorithm, described previously in the quasi-static electromagnetic simulation section. The motivation behind this simplification which allows an increase in speed of the calculation, allowing the increase the model's resolution and/or performing calculations for different magnetic field orientations in reasonable timescales. Allied with the CUDA implementation, this approach reduces the pressure calculation times by two orders of magnitude over that required for the full Navier-Stokes steady state solution.

The models
In this work a finite element geometric model is reported that approximates, to a first degree, the real geometry of the VS in order to answer the questions posed by the hypothesis of Roberts et al (2011). In this numerical simulation of the Lorentz force related pressure, some basic but biologically reasonable assumptions have to be made.
(1) The only currents present in the VS are ionic currents having their origin in the dark cells and flow towards the hair cells (Hibino and Kurachi 2006). (2) The ionic currents are present within the endolymph, with no leakage. Although marginal cells, transition cells and even epithelial cells may have a role in the maintenance of potassium balance (150 mM) in the endolymph, the current is assumed to come from specific dedicated dark cells located adjacent to the macula structures (Kimura 1969, Ciuman 2009).
(3) The current flow is local to the structures themselves and is balanced. This should be a reasonable assumption with the proviso that it is known that the saccule contains few if any dark cells and that a current may flow through to this chamber from the utricle but this is not modeled here. (4) The pressure is conserved inside the endolymph as there is no flexibility of the epithelium of the canals or utricle. This is a reasonable assumption for the normal working of the VS. (5) A true equilibrium state is assumed. Hence it is assumed that the time of exposure to a static magnetic field is longer than the 'long' time constant associated with cupula-canal operation . The 'long' time constants are of the order of a few seconds and hence the equilibrium state will be easily achieved. (6) The cupula deflection should only affect the transient values of the pressure until it reaches equilibrium, being a minor effect for the steady-state pressure calculation, so in this simulation the cupula is modeled as an infinite stiffness rigid plate (rather than a flexible flap). Additionally, estimates of the 'stiffness' term which relate pressure to deflection predict a deflection of 28 μm for pressures calculated with Kassemi et al's (2005) simulation. Oman et al (1987) predict an even lower value of only a few microns. Such values would correspond at most to one voxel in the resolution used in the simulation. Hence cupula deflections are not supposed to interfere with the steady state final calculation. (7) The cupula has the same conductivity as the endolymph. The cupula has a gel like constituency and is therefore considered as permeable to potassium ions as endolymph.
Hence current can flow unimpeded towards the hair cells but fluid cannot cross the boundary and a pressure difference is developed.
Two model geometries were constructed (figure 1), the first with a semicircular canal and an ampulla (CA), and the second similar, but adding the utricle (CAU). The first is designed to understand pressure propagation inside a simple geometry under ideal conditions, while the second model tries to approximate the problem to a more realistic situation.
Four current distributions are analyzed, two for each model, calculating the pressure difference on the cupula for each.
(1) Opposite current sources using the CA model ( figure 2(b)), parallel to the cupula, to produce a force perpendicular to the membrane in what would be a maximum theoretical limit. This situation reproduces Roberts et al's (2011) supplement simple calculation. (2) Ampullar currents, using the CA model, present in the surface of the crista ampullaris ( figure 1(a)).
(3) Utricular sources, using the CAU model, to measure the contribution of the long range current effects on the cupular membrane ( figure 1(b), excluding the ampullar sources and sinks). (4) Ampullar and utricular sources, using the CAU model, to measure the contribution of all current sources in the most realistic situation ( figure 1(b)). For all models and calculations, the magnetic field is oriented in the z direction unless stated otherwise. It is for this direction that a maximum pressure is expected, since the force components will be oriented orthogonally to the membrane.
The cupula and endolymph are considered to have the same conductivity of 1 S m -1 (i.e. the ions flow in the membrane as easily as in the endolymph), while the remaining tissue is non conductive (no ion fluxes through the tissues or into the perilymph). Since the whole volume has the same conductivity (σ ), this value has no influence on the calculation of the current densities. Total available ionic currents in the structures were considered as used by Roberts et al (2011) based on literature values, with a total of 0.7 μA for ampullar sources and 3.31 μA for the macula. A hair resting current of 100 pA was considered in preference to the most conservative figure of 10 pA used by Roberts et al. The higher figure appears to be entirely reasonable based on studies of single hair cells (Alexandrov et al 2001). Roberts et al (2011, supplement) concedes that a higher current than 10 pA per cell is flowing yet chooses the most conservative current in order to limit the calculated pressures to sensible values based on the simple model they used. Cochlea hair cells have even higher resting state currents than vestibular cells and these appear to be of the order of 1-4 nA (Johnson et al 2011).

Semicircular canal and ampulla
At this first stage it was prudent to model only one of the canals in order to understand the basic mechanism of pressure propagation without added complexity ( figure 1). An exact geometry would be preferable, but data on the exact orientation of all crista ampullaris and geometry of the endolymphatic ducts was not available. This first basic model was constructed based on elemental geometrical shapes that should represent the structure of a semicircular canal and ampulla. The ampulla is constituted by the crista, membrane, dark cells and hair cells. The dimensions we use are based on measures of the VS first performed by Curthoys and Oman (1987) and improved by Rother et al (2003). The semicircular canal is modeled by a torus with a minor radius of 0.2 mm and a major radius of 3.0 mm. The ampulla is considered to be an oblate spheroid of bigger radius 0.9 mm and small radius 0.7 mm. Inside the ampulla, the crista was assigned a Gaussian function with a height equal to 0.5 mm. The cupula was defined as a central region in the ampulla above the crista with 0.3 mm thickness and has an area of approximately 1.1 mm 2 . The size of the model is 323 × 338 × 73 voxels, with a resolution of 0.02 mm, and orientations are coincident with the field coordinates (z direction in the model corresponds to B z ).

Opposite currents.
Sources of ionic current (dark cells) were assigned to the bottom of the ampulla, close to the membrane, while the sinks (hair cells) were symmetrically assigned to the top side of the ampulla, with a total current of 0.7 μA flowing ( figure 2(b)). This value comes from measures of the number of hair cells (around 7000) and the resting current for each hair cell (100 pA). Having the current flowing parallel to the membrane should produce a Lorentz force perpendicularly oriented to the membrane, in what would be a theoretical maximum as calculated by Roberts et al (2011).

Crista currents.
Sources of ionic current were assigned to the bottom of the crista, while the sinks were assigned to the top of the crista, mimicking the behavior of hair cells in ionic flow (figure 2(a)). In both scenarios a total current of 0.7 μA is produced from the dark cells and flows into the hair cells.

Semicircular canal, ampulla and utricle
The second model was designed to test the effect of the ionic current originated in the utricular macula on the cupula, using an ellipsoid as shown in figure 1. The dimension of the semiprincipal axes for the ellipsoid are 1.36 mm for x, 1.55 mm for y and 0.67 mm for z (Rother et al 2003). Total volume for the utricle is 5.9 mm 3 . The size of the model is 353 × 338 × 73 voxels, at a 0.02 mm resolution.

Utricular currents.
The utricular macula was modeled by setting the bottom of the utricle as currents sinks, with a surrounding ring of current sources, as shown in figure 1(b). The total current flux is 3.31 μA (corresponding to 33 100 hair cells, with a resting current of 100 pA each). It was designed to measure the effect of long range utricle currents on the cupula via the coupling of the pressures through the endolymph and canals and is assumed to be L-shaped, on the outer surface of the utricle (Curthoys et al 2009).

Utricular and ampullar sources.
This current configuration is created by simultaneously modeling the currents of the crista and the utricle. Utricle currents are not expected to have a peak pressure contribution for a z oriented field due to positioning in a different plane.

Results
The simplified pressure model simulations were compared to the solution of the full Navier-Stokes equations using the software package Nast3DGP, developed by the research group in the Division of Scientific Computing and Numerical Simulation at the University of Bonn (Griebel et al 1998), with a small modification to allow the introduction of force density as a function of position (f). We used a Reynolds number of 1000, QUICK differencing scheme, second order Adams-Bashfort time evolution, BiCGStab solver with 1000 maximum iterations using a minimum residual of 10 −6 . Current densities necessary to pass the Lorentz force to Nast3DGP were calculated using the quasi-static simulation mentioned above. The N-S simulation continued until a steady-state solution was found. Pressure differences across the membrane differed by no more than 15% compared to the approximated equations, fluid velocities were found to be no more than 3 × 10 −4 ms −1 , and the convective acceleration v · ∇v, was always less than 2 × 10 −10 ms −2 corroborating the validity of the approximation. The N-S simulation took 8 h to run on an Intel i7 950 3.06 GHz, in comparison with 10 min for the approximated equations on the same CPU and 1 min for the GPU version on a nVidia Geforce GTX 560 Ti card. The trans-cupular pressure, which is the difference between the average pressure of the fluid voxels adjacent to each side of the membrane, is calculated. A negative value of the pressure means that the pressure is lower at the right side of the membrane. Changing the thickness of the membrane between 0.1 and 0.3 mm does not influence the pressure calculation by more than 1%. This small variation is due to small changes in the force density field close to the cupula.

B z field
For the first model (CA), with sources and sinks placed to maximize the force exerted on the cupula, the trans-cupular pressure was calculated to be −3.3 mPa ( figure 3(a)). This is not a realistic physical model, and is only an indication of the maximum pressure for the given ionic current strength. Ampullar currents with both source and sink in the crista (figure 3(b)) are significantly smaller (−0.9 mPa), since they do not follow a parallel path with the cupula. For the second model (CAU), with utricular currents (figure 3(c)), the pressure difference came to −0.7 mPa, and when both utricular and ampullar currents were considered (figure 3(d)), the pressure was −1.6 mPa.

B x and B y fields
The same simulation was re-run for magnetic fields oriented in the x and y directions. The pressure differences for the four situations are written in table 1.

Discussion
It was possible using finite element models for electromagnetic and fluid equations to calculate a more realistic value of the trans-cupular pressure in the VS due to a Lorenz force field. The difference in results is significant and justifies the use of such detailed approach. Even though the reduction is considerable, it is still enough to justify the MFIV and nystagmus effects reported in the vicinities of large magnetic fields. The summary values in table 1 show that the pressure on the membrane is strongly dependent on the orientation, as the force field moves the regions of high pressure from the membrane to the walls of the VS. Also, the existence of components in the force field with opposite directions leads to its cancellation, causing a large drop in the trans-cupular pressure. Even in the biologically unrealistic scenario of opposite current sources producing a quasiparallel current (figure 3(a)), there is still a considerable amount of cancellation which justifies the magnitude in pressure of 3.3 mPa, in comparison with the 5 mPa calculated by Roberts et al (2011). The pressure drop due to these cancellations is more than one order of magnitude for the more realistic geometries.
For the CAU model, the utricular currents (3.31 μA) are significantly larger than the available ampullar currents (0.7 μA), but due to the macula's 'L-shape' (Curthoys et al 2009), there are force components in all directions. This causes the utricle originated pressure (∼0.7 mPa) to be largely dissipated, but present in all orientations, with a contribution of at least 45% of the total pressure. Not surprisingly, the maximum pressure derived from the utricular sources does not necessarily occur for a pure B z magnetic field, since utricular current flow is a function of the surface of the macula, which in turn is not parallel to the cupula.
The pressure for the combined utricular and ampullar currents is verified to be the sum of the individual sources, reaching a maximum value of −1.6 mPa for the z-orientation.
The pressure results will scale linearly with the field strength, since B only intervenes as a constant during the Lorenz force calculation. For example, at 3 T, pressures would be 3/7 of the presented values. This linear effect is evident in the experimental data shown in Roberts et al (2011).
Estimates of trans-cupular pressure causing nystagmus responses are available from simulations of caloric stimulation and provide useful comparisons. Kassemi et al (2005) consider a 2D finite element calculation of caloric induced trans-cupular pressure. They calculate that a temperature increase of 1 • C induces a pressure differential of ∼4 mPa, not far away from the 3.3 mPa estimated by Valli et al (2003). The measured slow-phase nystagmus due to the standard caloric tests, which induce a temperature change in the lateral semi-circular canal of approximately 1 • C (Cawthorne and Cobb 1954), is approximately 20 • s -1 either with eyes open in darkness (Munro and Higsonaf 1996) or with eyes closed (Coats and Smith 1967). Roberts et al (2011) measure a typical slow phase velocity of 10 • s -1 inside the 7 T magnetic field. Assuming the same mechanism (transcupular pressure) is causing the nystagmus and a linear response between trans-cupular pressure and slow phase velocity, a figure of 1.6-2.0 mPa due to a 7 T field may be expected by combining the above observations, in accordance with the assumptions and results presented in the current work.

Conclusions
Using a simple but realistic model of the inner ear it is possible to calculate the pressures on the cupula due to Lorentz forces originated by ionic currents in a static magnetic field. The values of the pressures are substantially lower than the ones calculated by assuming a purely perpendicular component of the force against the cupular membrane. Even though the pressure is dissipated throughout the fluid, there is still an appreciable trans-cupular pressure capable of justifying the nystagmus in a 7 T magnetic field. Using this more realistic model supports the hypothesis that the Lorentz forces due to endolymph currents and applied magnetic fields may be used to explain the phenomena of magnetic field induced vertigo. Consistent with available literature, a more reasonable hair-cell resting current of 100 pA per unit has been assumed. The most influential factor for the trans-cupular pressure is the component of the force field orthogonal to the membrane, meaning that field direction, geometry and localization of ionic sources and sinks are the dominant factors to be taken into account. Accurate modeling of the VS is crucial to calculate the effect of the current distributions on the cupula. Better knowledge of vestibular current geometries (including better measurements of their magnitude) would be vital in order to refine the model described here. Better knowledge of utricle and saccule geometry and current flow would also increase the validity of the model. However, even with current understanding it should be possible to model a full head system with both the inner-ear vestibular systems and hence all 6-ampullae. In this way the exact response of the human to magnetic field strength and orientation can be modeled and compared to human perception and response.
Dropping the convective acceleration term in the fluid equations is shown to be acceptable when the force densities are large enough. This approximation is useful for simulating models with higher resolution, or for calculating multiple Lorentz force orientations quickly and easily.