Non-invasive imaging of neural activity with magnetic detection electrical impedance tomography (MDEIT): a modelling study

Objectives. (1) Develop a computational pipeline for three-dimensional fast neural magnetic detection electrical impedance tomography (MDEIT), (2) determine whether constant current or constant voltage is preferable for MDEIT, (3) perform reconstructions of simulated neural activity in a human head model with realistic noise and compare MDEIT to EIT and (4) perform a two-dimensional study in a saline tank for MDEIT with optically pumped magnetometers (OPMs) and compare reconstruction algorithms. Approach. Forward modelling and image reconstruction were performed with a realistic model of a human head in three dimensions and at three noise levels for four perturbations representing neural activity. Images were compared using the error in the position and size of the reconstructed perturbations. Two-dimensional MDEIT was performed in a saline tank with a resistive perturbation and one OPM. Six reconstruction algorithms were compared using the error in the position and size of the reconstructed perturbations. Main results. A computational pipeline was developed in COMSOL Multiphysics, reducing the Jacobian calculation time from months to days. MDEIT reconstructed images with a lower reconstruction error than EIT with a mean difference of 7.0%, 5.5% and 11% for three noise cases representing current noise, reduced current source noise and reduced current source and magnetometer noise. A rank analysis concluded that the MDEIT Jacobian was less rank-deficient than the EIT Jacobian. Reconstructions of a phantom in a saline tank had a best reconstruction error of 13%, achieved using 0th-order Tikhonov regularisation with simulated noise-based correction. Significance. This study demonstrated that three-dimensional MDEIT for neural imaging is feasible and that MDEIT reconstructed superior images to EIT, which can be explained by the lesser rank deficiency of the MDEIT Jacobian. Reconstructions of a perturbation in a saline tank demonstrated a proof of principle for two-dimensional MDEIT with OPMs and identified the best reconstruction algorithm.

1. Introduction 1.1.Background Magnetic detection electrical impedance tomography (MDEIT) is a novel non-invasive imaging technique built upon the principles of electrical impedance tomography (EIT) and magnetometry (Chen et al 2020).The working principle of EIT is to attach an array of electrodes to the boundary of a region of interest, inject an alternating current (AC) through pairs of electrodes and measure the voltage on all non-injecting electrodes.This is done before and during a local, internal change in electrical impedance and the difference in voltage between the time points is used to reconstruct an image of the local impedance change (Adler and Holder 2021).EIT can be used to image functional neural activity because there is a local change in impedance of ;1% associated with the neuronal depolarization of an action potential (Cole and Curtis 1939, Holder 1992, Liston et al 2000, Liston 2003, Tarotin et al 2019).
In the nerve, EIT has successfully been used to image the fast neural activity in the sciatic and vagus nerves of the rat and pig respectively, allowing for the fascicular organisation of the pig vagus to be determined for the first time (Aristovich et al 2018, Ravagli et al 2020, Thompson et al 2022).In the brain, EIT has been deemed impractical for neural imaging with scalp electrodes as the signal-to-noise ratio (SNR) was too low (Holder 1989, Gilad andHolder 2009).However, EIT has achieved a resolution of 200 μm and 2 ms in the rat brain with an array of epicortical electrodes, which was limited to the cortex (Aristovich et al 2014, 2016, Faulkner et al 2018b).
The limitations of fast neural EIT with scalp electrodes are largely attributable to the skull, which is more electrically resistive than the surrounding tissue and can attenuate the signal by a factor of 10-100 (Liston 2003, Romsauerova et al 2006, Gilad et al 2015).In order to overcome this, AC can be injected with scalp electrodes as in EIT and the magnetic field outside around the head can be measured instead of the voltage on the scalp.Whilst the injected current is still attenuated by the skull and shunted by the scalp, the magnetic field is not attenuated by the skull (Singh 2014), meaning magnetic measurement could theoretically produce an increase in the quality of the reconstructed image.
MDEIT was first demonstrated in a saline tank using a rubber cylinder as a perturbation.AC was injected at 16 Hz through two electrodes and the magnetic field was measured at 12 positions using superconducting quantum interference (SQUID) magnetometers.Images of the current distribution were reconstructed and a minimum-norm estimate was used to localise the position of the perturbation (Ahlfors and Ilmoniemi 1992).Magnetic search coils at 240 locations have also been used in conjunction with 10 kHz, 100 kHz and 1 Mhz AC at 10 mA through two electrodes in a saline tank to reconstruct images of a resistive object, which was then extended to images of lung inspiration and expiration in a human.The reconstructed images qualitatively corresponded to the perturbation; however, a quantitative spatial resolution was not stated in either case (Tozer et al 1999, Ireland et al 2004).Forward modelling of fast neural MDEIT has previously been performed for activity in the human brain (Ahadzi et al 2004, Gilad et al 2009), which concluded that the SNR is comparable to that of EIT; this was followed by measurements of the MDEIT signal with scalp electrodes in humans which agreed with the results of the forward modelling but indicated that an experimental time of three hours would be required for image reconstruction due to low SNR (Gilad 2007).For this reason, image reconstruction was not performed (Ahadzi et al 2004, Gilad et al 2009).Direct current (DC) was used in this study; however, it is now known that the EIT SNR is larger at a higher frequency, peaking at 1475 Hz for neural activity in the brain (Faulkner et al 2018a).It is not known whether the SNR for MDEIT will be largest at the same frequency, but it can be inferred that a larger SNR can be expected than was measured at DC Overall, there is a limited body of literature on MDEIT, with approximately 15 publications, most of which only consider two-dimensional imaging and none of which perform three-dimensional MDEIT for neural imaging.
Since these studies, there have been significant advances in the optimisation of finite element models (FEMs), AC injection protocols and the AC injection frequency for EIT that have direct applications in MDEIT (Aristovich et al 2016, Jehl et al 2016, Faulkner et al 2017, 2018a).There has also been progress in magnetometry, with new commercial and research-level SQUID magnetometers increasing in sensitivity (Neuromag 2008, Fedele et al 2015, Storm et al 2017, CTF 2021) to ∼1 fTHz −1/2 and decreasing in cost and operating temperature by cooling the systems with liquid nitrogen instead of liquid helium (Faley et al 2017).Optically pumped magnetometers (OPMs) have now been developed with a sensitivity of ∼10 fTHz −1/2 , whilst operating at room temperature (Shah and Wakai 2013, Tierney et al 2019, Quspin 2022).The benefit of OPMs over SQUID magnetometers is that they are less expensive to operate, housed in a small form factor and are portable, whereas SQUIDs are usually housed in a large dewar and are fixed in place, this has led to a significant uptake of OPMs for MEG measurements (Hill et al 2020, Seymour et al 2021).

MDEIT forward and inverse problems
The forward problem of MDEIT is to calculate the magnetic field at an arbitrary position outside a conductive medium given the current on the injection electrodes and the electrical conductivity distribution inside the body.The computational implementation is usually performed using FEMs, which coarse-grain the problem by discretising the space into voxels across which the conductivity is constant (Bathe 2007).A quasi-static approximation can be used since the frequencies under consideration are generally in the kHz range which is below the ∼1 MHz limit above which the approximation is no longer valid (Zhang and D Li 2014).
To solve the MDEIT forward problem, the EIT forward is first solved to find the voltage and current distribution in the region of interest (Adler and Holder 2021), once this is known, the magnetic field at any point  r can be calculated using the Biot-Savart law for a FEM where   B r ( ) is the magnetic field at  r , μ 0 is the magnetic permeability of free space,  J i is the current density in the ith element,  ¢ r i is a vector from the centre of the ith element to the magnetic sensor and V i is the volume of the ith element (Jackson 1999).Equation (1) shows that the magnetic field follows an inverse square law, meaning that the size of the MDEIT signal will decrease with the square of the distance from the source (i.e. the perturbation) to the sensor.
The forward problem can be expressed in matrix form as is the magnetic field at the magnetometers and  is the forward operator where N m is the number of measurements and N e the number of elements in the FEM.For small changes in conductivity δσ and magnetic field δb, the forward operator can be linearised as where J is the Jacobian matrix (or sensitivity matrix) relating the change in the magnetic field at the ith sensor to the change in conductivity of the jth element in the FEM (Adler and Guardo 1996).Since the magnetic field is a vector quantity, MDEIT can be performed with the measurement of one, two or three components of the magnetic field (figure 1) and the Jacobin's for each component can be concatenated to form one Jacobian for all measurements.The Jacobian can be calculated using the adjoint state method in COMSOL Multiphysics (COMSOL Multiphysics 2015).
The inverse problem of MDEIT is to calculate the distribution of conductivity changes in the conductive region given the magnetic field change b at the magnetometers.The problem is nonlinear, ill-conditioned and ill-posed so the problem is linearised and regularised in order to find a solution (Lionheart 2004, Hansen 2010).Linearisation is done by only considering small changes in conductivity and regularisation is done by incorporating a priori information into the solution (Holder and Khan 1994, Adler and Guardo 1996, Lionheart 2004).Since J is, in general, a rectangular matrix with more columns than rows, the inverse of J is calculated using the singular value decomposition and the Moore-Penrose generalised inverse (Penrose 1955, Hansen 2010), the solution can then be expressed as where  is the regularisation matrix, λ is the regularisation parameter, and σ λ is the reconstructed conductivity distribution.For the case of Tikhonov regularisation L = I, the identity matrix (Phillips andTechnique 1962, Tikhonov 1963).For NOSER regularisation, L = diag(J T J) (Cheney et al 1990).The optimal value of λ can be calculated using methods such as heuristic selection (Graham and Adler 2006), the L-curve method (Hansen 1992(Hansen , 1998)), the fixed noise figure (Adler and Guardo 1996), the BestRes method (Graham and Adler 2006) and generalised cross-validation (GCV) (Hansen 1998.GCV is a method that seeks to match σ λ with the measured data b as well as possible and has been dubbed 'the statistician's choice' (Hansen 2010).This is done by minimising the GCV function 1.3.Noise-based correction Noise-based correction (NBC) is a post-processing method for suppressing the output of voxels in the reconstruction as a function of their contribution to the measured noise, with a larger noise contribution corresponding to a greater suppression (Aristovich et al 2014, Faulkner et where i represents the ith component of the quantity and i ä {1, 2, K, N e }.

Purpose
The main contributions of this work are the development of a computational scheme for three-dimensional MDEIT in anatomically realistic domains, the first performance of three-dimensional MDEIT for fast neural imaging and comparison with EIT, the first demonstration of MDEIT with an OPM and the first comparison of Tikhonov regularisation with NOSER (with and without NBC) for MDEIT.The purpose of this work was to answer the following questions: 1. Can an efficient computational pipeline for MDEIT be implemented?
2. Is constant current or constant voltage injection preferable for MDEIT?
4. Does MDEIT work in a saline tank with an OPM and what is the best reconstruction algorithm?

Experimental design
This study can be separated into three branches: the development and implementation of a computational scheme for MDEIT, the forward and inverse modelling of MDEIT and EIT and the demonstration of MDEIT in a saline tank with an OPM.

Algorithms
For the calculation of the MDEIT Jacobian, the forward method is computationally impractical for large meshes (millions of elements).Therefore it is imperative that the adjoint state method be implemented for the calculation.The adjoint state method is regularly used in EIT (Polydorides and Lionheart 2003), magnetic induction tomography (MIT) (Soleimani and Lionheart 2006) and in magnetic and electric impedance imaging for geophysical applications (Chen et al 2005, Plessix 2006, Dorn et al 2008); however, it has never, to the best of the authors' knowledge, been explicitly implemented for the case of MDEIT in the quasi-static frequency range.

Computational model
The FEM used in this work was an anatomically accurate representation of a human head comprising seven different tissue types.This was essential for the comparison of MDEIT with EIT since it is the attenuation of the signal by the skull that was considered to be the primary factor that makes fast neural EIT infeasible with scalp electrodes (Ahadzi et al 2004, Gilad andHolder 2009).The perturbations used were designed to accurately represent neural activity in the brain, considering a local impedance change of 1% (Holder 1992, Liston et al 2000, Liston 2003, Tarotin et al 2019) and a volume of 3.86 cm 3 (Pastor et al 2003) at four different locations approximately corresponding to the cortex, cingulate gyrus, thalamus and pons (Nowinski 2011).This is consistent with perturbations considered in fast neural EIT (table1) (Aristovich et al 2016).
1.5.3.Software COMSOL Multiphysics (COMSOL AB 2022) was the software chosen for the calculation of the MDEIT Jacobian.This was due to the software's commercial availability, ease of use and built-in 'sensitivity' functionality which allows for the calculation of the sensitivity of the magnetic field with respect to the conductivity.To calculate the sensitivity, COMSOL requires that a geometry and material be defined in the region surrounding the conductive region, so a cubic region of air surrounding the head was incorporated into the model.
1.5.4.Current/voltage injection Constant current and voltage injection were compared in order to assess whether constant voltage injection would produce an increase in the SNR over constant current injection, which is standard in EIT.When constant current injection is performed, the distribution of current density inside the volume can change with an impedance perturbation, but the total current remains the same.Since the magnetic field is proportional to the current density (equation ( 1)), constant voltage injection could increase the SNR, since it would allow the total current to change as well as the distribution.

Injection protocol
The injection protocol used was selected as the protocol that maximised the current density in the brain (grey and white matter voxels) (Faulkner et al 2017), this region was selected instead of the individual perturbations since in practice, it is likely that the location of the neural activity will not be known.

Magnetometer sensitivity
The magnetometers considered in the computational modelling were not considered to correspond to any physical magnetometer in particular but served to represent a typical magnetometer which would be available in practice.(2) The expected resolution can be approximated by calculating the size of voxels in the hexahedral mesh such that the number of voxels equals the number of measurements (i.e. one measurement for each voxel value).Since perturbations of 20 mm diameter were considered, a resolution of at least 20 mm was necessary.The resolution can therefore be approximated as where D is the mean dimension of the model and N m is the number of independent measurements.For a resolution of 20 mm N m must be 903 giving the number of measurement electrodes/magnetometers as 30 (for 31 injection electrode pairs).However, Jacobian matrices are often rank-deficient in impedance imaging, effectively reducing the number of independent measurements (Luppi Silva et al 2017).For this reason, it makes sense to exceed the minimum number of measurements to ensure this condition is met.64 is a suitable number of magnetometers/measurement electrodes which can be expected to comfortably meet this condition (and condition (1)), taking into account rank deficiency (Luppi Silva et al 2017).(3) The computational resources required to compute the forward solution and Jacobian increase as the number of magnetometers and electrodes increases, therefore this number must be chosen such that conditions (1) and (2) are met without exceeding the computational capabilities available.64 was the largest number for which this was possible.(4) The practicalities of experimentation must be taken into account when choosing the number of electrodes and magnetometers, QuSpin OPMs (Quspin 2022) and/or Cerca OPM systems (Cerca 2023) are technologies with future potential to enable MDEIT with OPMs, both of which have commercial on-scalp systems for 64 magnetometers.It is worth noting that any number of electrodes and magnetometers can be chosen in principle and a full optimisation study can be performed to accurately determine the best number of electrodes/magnetometers; however, this was not deemed necessary for the purposes of this work.The relationship between the number of electrodes was chosen such that the number of measurement locations was the same for EIT and MDEIT.
1.5.8.Noise Three noise cases were considered in this work to capture the expected SNR and reconstructed image quality associated with different, realistic, hardware.The noise was split into the multiplicative (i.e.current source noise) and additive (i.e.environmental/magnetometer noise) components.The electric noise for the first noise case corresponded to the measured noise with electrodes on the scalp of a human (supplementary material).The magnetic noise was calculated from literature values of the environmental noise in a magnetically shielded room (Storm et al 2017), the sensor noise of a modern OPM or SQUID system (Neuromag 2008, CTF 2021, MAG4Health 2021, Quspin 2022) and the measured current source noise.The second noise case considered the same additive noise with a 56-fold reduction in current source noise, this was considered to be a realistic improvement that is achievable in the next five years of hardware development using CMOS-based architectures (ROHM SEMICONDUCTOR 2023) in combination with advanced denoising techniques (Dos Reis Filho 2022).The final noise case considered the same multiplicative noise as noise case 2 with a 1000-fold reduction in the magnetometer noise.This is a potential, although ambitious, future improvement in magnetometry approximately corresponding to the limit of OPM sensitivity (Savukov et al 2005) or a 10-fold decrease in the noise of the most sensitive SQUID magnetometers (table 3)  1.5.9.Tank study For the tank study, one OPM was used and sequentially placed in 25 locations around the tank to simulate an array (Chen et al 2020).This was a practical limitation as only one OPM was available for the study and was not a design choice.Due to this limitation, the tank study cannot serve as a full validation of fast neural MDEIT but serves as a proof of principle for MDEIT with an OPM and to compare reconstruction algorithms using measured data.The injection frequency was chosen as 90 Hz because the OPM used was tuned to frequencies in the 3-100 Hz range, and it was expected that noise would be lower at the higher end of the frequency range (Faulkner et al 2018a, Quspin 2022).The reconstruction algorithms considered for the tank study were chosen based on their applicability to fast neural imaging.Only one-step regularisation methods were considered because they are faster to implement than methods such as total variation regularisation (Hao et al 2014).The effect of NBC was studied because it has never been applied to MDEIT before.

Computational model
An anatomically realistic 3D FEM comprising 3.2 M tetrahedral elements and seven different regions of electrical conductivity, each representing a different tissue, was used for all simulations (table 2).This was constructed and segmented from MRI and CT images of a real human head (Jehl et al 2016).For the Jacobian calculation, the same FEM was considered inside a cubic region of air, the combined model comprised 4.4 M tetrahedral elements.For Image reconstruction, a hexahedral mesh comprising 375k cubic elements was used.
32 scalp electrodes of 5 mm radius and 1 kΩ contact impedance were placed in the EEG 10-20 standard positions on the scalp of the FEM and an additional 34 electrodes of 5 mm radius were placed on the scalp in an approximately symmetric configuration by eye (Jasper 1958).All electrodes were used for voltage measurement, but only the first 32 were used for AC injection.A ground node was placed at the nape of FEM and was used as a reference for voltage measurements.64 magnetometers were considered to be in a helmet shape above the scalp of the FEM at a distance of 7 mm from the surface of the scalp and all three components of the magnetic field were calculated at each magnetometer location.
The AC injection protocol was selected such that the current density was maximised in the region of interest, which was the whole brain (Faulkner et al 2017). 1 mA was maintained on the injection electrodes for the case of constant current injection and the voltage which maintained a current of 1 mA on injection electrodes for the unperturbed case was applied to the injection electrodes for constant voltage injection.
Four approximately spherical perturbations of 1% local increase in conductivity were considered at four depths in the brain, the volume of each perturbation was 3.86 cm 3 and only included elements of the FEM that corresponded to white or grey matter (table 1 and figure 2) (Liston 2003, Aristovich et al 2016).
The FEM computational implementation for MDEIT was validated by comparing the solution for the magnetic field with the semi-analytical solution of current flowing through a long wire (Charitat and Graner 2003).A FEM comprising 138k elements in the geometry of a wire of length 200 mm and diameter 0.5 mm was meshed with one electrode at either end.Currents of 1 A, 2 A and 3 A were simulated to be flowing through the wire and the magnetic field was calculated at 100 radial positions from the wire ranging from 1 to 10 mm from the centre of the wire.

Noise
Three noise cases were considered for the forward and inverse modelling, distinguishing between the additive and multiplicative components (table 3).These corresponded to the current state of measured/calculated noise (noise case 1), reduced current source noise from 0.058% to 0.001% (noise case 2) and the same reduced current source noise combined with reduced magnetometer noise from 10 to 0.01 fTHz −1/2 .The noise figures in each case were scaled according to the number of measurement averages deemed realistic in one hour of recording, for 31 injection pairs and an evoked activity duration of 500 ms, this was 232 averages (table 3).For noise case 1, the electric noise was considered equal to the measured additive noise and multiplicative noise (supplementary material).The multiplicative magnetic noise was taken to be equal to that of electric measurements and the additive magnetic noise was calculated by estimation of the intrinsic sensor noise and  environmental magnetic noise at ∼1.5 kHz in a magnetically shielded room.The intrinsic sensor noise was taken to be 10 fTHz −1/2 , corresponding to the noise of a typical SQUID magnetometer system or modern OPM (CTF 2021, Quspin 2022) and the environmental noise was taken to be 0.2 fTHz −1/2 as was measured in a magnetically shielded room (Storm et al 2017).

Algorithms
For all forward solutions, the interior current distribution and electrode voltages were calculated using the Electrical Impedance and Diffuse Optical Reconstruction Software (EIDORS) in Matlab and the magnetic field was calculated using custom-written code in Matlab (Polydorides and Lionheart 2002, MATLAB 2021, Polydorides et al 2022).The MDEIT Jacobian was calculated with a sensitivity study in COMSOL Multiphysics utilising the adjoint state method and calculation of the EIT Jacobian was performed using the calc_jacobian function in EIDORS (COMSOL AB 2022, Polydorides 2022).All simulations were performed in three dimensions.
For image reconstruction, 0th-order Tikhonov regularisation (TR) with simulated NBC was used to reconstruct all images in three dimensions and leave-one-out GCV was used to find the regularisation parameter.Forward simulations were performed for EIT and 3-axis MDEIT and reconstruction was performed for EIT, 1-axis MDEIT and 3-axis MDEIT.Additive and multiplicative noise were added to the simulated signal before being fed to the reconstruction algorithm (table 3).

Mesh convergence
Six FEMs representing the 3D head model (each comprising seven tissue types) (figure 2), each with a different number of elements, were meshed using COMSOL Multiphysics (COMSOL AB 2022).For each mesh size, three FEMs were meshed of that size, varying by a comparatively small number of elements.An MDEIT forward solution was computed for each mesh and the mean of the largest 10% of the magnetic field was considered as a function of the number of elements in the mesh.The convergence error was calculated as the mean difference in the magnetic field between mesh n m and n m+1 where n m ä {1, 2, 3, 4, 5, 6}, which gave five mesh refinement steps.If the mesh converged on mesh refinement step n s , then mesh number n m = n s was taken to be the coarsest mesh for which convergence had been achieved.The mesh variability was calculated as the mean difference in the magnetic field between all permutations of the three meshes at each refinement step.Mesh convergence was assessed by heuristic inspection of the mesh variability and convergence error.

Data analysis
The SNR of the largest single raw change and the mean of the largest 10% of raw changes were calculated from the forward modelling and used to calculate the SNR which was defined as where the change is either a magnetic field or voltage change.
The images were compared using the total reconstruction error of the image, defined as where E p was the position error and E V was the volume error of the reconstruction.The position error was defined as where  x true was the true location of the centre of mass of the perturbation,  x recon was the location of the centre of mass of the reconstructed perturbation of the thresholded image and x FEM was the mean dimension of the FEM.The volume error was defined as where V true was the true volume of the perturbation, V recon was the volume of the reconstructed perturbation of the thresholded image and V FEM was the total volume of the FEM.The error was analysed using repeated measures ANOVA and multiple comparison tests to assess the significance of the difference in total reconstruction error between cases.100 reconstructions were performed for each case.

Tank study
A tank of 80 mm in diameter and 70 mm in height was 3D printed with 16 equally spaced recesses on the interior wall at a height of 35 mm for Ag/AgCl electrodes of 9 mm in diameter (Formlabs 2022).One Quspin QZFM Gen 2.0 OPM was used to measure the vertical component of the magnetic field only (Quspin 2022).The OPM was held at a radius of 71 mm from the centre of the tank in the same plane as the electrodes and was controlled by a system of two gears and a stepper motor, allowing the OPM to be rotated in a plane around the tank.A plastic cylinder of 25 mm in diameter was used as a conductivity perturbation and was placed 20 mm from the centre of the tank.The tank was placed inside a 3-layer magnetically shielded chamber comprising 2 layers of Mumetal®and one layer of aluminium (figure 3).
The experimental procedure was executed in the following steps: (i) Start with the tank containing only saline.
(ii) Position the OPM in the first location.
(iii) Perform the entire injection protocol.
(iv) Move the OPM to the next position and repeat step (iii).
(v) Repeat step (iv) for all OPM positions.
AC was injected in a 'skip 2' protocol (i.e. between electrodes [1, 4], [2, 5], K, [14, 1]) at 90 Hz for 1s per injection pair.The current was injected at a peak-to-peak amplitude of either 0.264 mA, 0.8 mA or 2.4 mA which were paired with amplifier gains of 3×, 1× or 0.33× respectively.The magnetic field data was output as a voltage which was sampled by an analogue to digital converter at a sampling frequency of 5 kHz.17 datasets were collected in total.The raw data was filtered with a 3rd order Butterworth bandpass filter with a bandwidth of ±5 Hz, centred at 90 Hz.The data was then demodulated using the Hilbert Transform and filtered with a 3rd order Butterworth lowpass filter with a cut-off frequency of 5 Hz.A magnetic field change between the perturbed and unperturbed cases was calculated by taking the difference of the mean of the magnetic field over the middle 20% of the injection time.The MDEIT Jacobian was calculated using the forward method and images were reconstructed on a 2D FEM comprising 2014 triangular elements using six different reconstruction algorithms: 0th-order Tikhonov regularisation (TR), 0th-order TR with simulated noise-based correction (NBC), 0th-order TR with real NBC, NOSER, NOSER with simulated NBC and NOSER with real NBC (Cheney et al 1990).The regularisation parameter was found using leave-one-out cross-validation for each algorithm.The reconstructed images were thresholded at 50% the largest negative change in conductivity and the total reconstruction error was calculated as in (9), with volume error becoming area error in 2D.

Mesh convergence and model validation
The model validation resulted in a mean difference of 0.17% between the analytical solution and numerical solution across all current levels and measurement positions.There was a linear correlation with R 2 = 1.00 between the analytical and numerical solutions.
The mesh convergence was deemed to have been achieved at mesh refinement step four (figure 4).

Constant current versus constant voltage injection
For the comparison of the SNR between constant current and voltage injection, only noise case 1 was considered with 3-axis MDEIT.The SNR of the largest change was larger for constant current injection for perturbations 1, 2 and 4 with a mean decrease of 11% when using constant voltage injection.The SNR of the mean of the largest 10% of changes was larger for all perturbations for constant current injection, with a mean decrease of 0.79% for constant voltage injection (figure 5).

EIT versus MDEIT
Analysis of the regularisation parameters showed that the regularisation parameter used was larger for EIT than 1-axis and 3-axis MDEIT for an equivalent SNR (figure 6).On visual inspection, image reconstructions of the conductivity perturbation showed a correspondence between the reconstructed perturbation's size and location and that of the true perturbation (figures 1 and 7), the image quality decreased as the noise increased and as the perturbation depth increased (figure 7).On visual inspection, it was also concluded that EIT produced images of an inferior quality to MDEIT (1-axis and 3-axis), with MDEIT images having fewer artefacts and a clearer boundary between the reconstructed perturbation and the rest of the FEM (figure 7).
For noise case 1, MDEIT had a smaller SNR than EIT when considering the mean of the largest 10% of changes and a larger SNR than EIT for perturbations 1, 2 and 3 when the single largest change was considered (figure 8(a)).For noise case 2, MDEIT had a smaller SNR than EIT for all perturbations when considering either the largest change or the mean of the largest 10% of changes (figure 8(b)).For noise case 3, MDEIT had a larger SNR than EIT for all perturbations when considering either the largest change or the mean of the largest 10% of changes (figure 8(c)).The SNR of the MDEIT signal decreased more rapidly as a function of the depth of the perturbation than the SNR of EIT.
For noise cases 1, 2 and 3, EIT reconstructed images with a significantly larger total reconstruction error than 3-axis MDEIT for all four perturbations (P < 0.001, multiple comparison test, N = 100).EIT reconstructed images with a larger total reconstruction error than 1-axis MDEIT for all perturbations and noise cases which was significant for all perturbations for noise case 3 (P < 0.001, multiple comparison test, N = 100) and perturbations 1, 2 and 4 for noise cases 1 and 2 (P < 0.05, multiple comparison test, N = 100) (figure 9).
A rank analysis of the Jacobians for EIT, 1-axis MDEIT and 3-axis MDEIT showed that all three techniques had rank-deficient matrices and the degree of rank deficiency was larger for EIT than 1-axis MDEIT and 3-axis MDEIT (table 4).

MDEIT in a saline tank
The regularisation parameter used was unique to each reconstructed image but was consistently lower for TR than NOSER.Since NBC is a post-processing technique, it had no effect on the regularisation parameter (figure 10).For all reconstruction algorithms, the location and size of the reconstructed perturbation could be seen to correlate with the true location and size of the perturbation (figure 11).TR with simulated NBC had a lower total reconstruction error than all other reconstruction algorithms with a value (mean ± SE) of 12.2% ± 2.5%, comprising 9.43% position error and 2.72% area error.This was significant with respect to TR, and NOSER with real NBC (P < 0.05, multiple comparison test, N = 17) and insignificant with respect to TR with real NBC (P = 0.12, multiple comparisons test, N = 17), NOSER (P = 0.95, multiple comparison test, N = 17) and NOSER with simulated NBC (P = 1.0, multiple comparison test, N = 17) (figure 12).

Summary of results
A FEM computational scheme for the MDEIT forward problem was successfully developed and implemented for the case of an anatomically realistic human head and was verified with respect to an analytical solution.In addition to this, an efficient method for the calculation of the MDEIT Jacobian was introduced and implemented, utilising the adjoint state method in COMSOL Multiphysics (COMSOL AB 2022).
From the forward modelling, it was concluded that constant current is superior to constant voltage injection, the SNRs of EIT and MDEIT are similar given current noise, EIT's SNR is larger if the current source noise is reduced 56-fold and MDEIT's SNR is larger if the magnetometer sensitivity is increased 1000-fold.
Reconstructions of four regions of simulated neural activity perturbations across three noise cases showed that 1-axis MDEIT consistently reconstructed superior images to EIT with a mean difference in reconstruction error of 7.0%, 5.5% and 11.0% across all perturbations for each noise case respectively.3-axis MDEIT reconstructed superior images to 1-axis MDEIT in 11 of 12 perturbation and noise combinations with a mean difference in reconstruction error of 4.8%, 1.6% and 0.71% across all perturbations for each noise case respectively.The EIT Jacobian was calculated to be more rank-deficient than the 1-axis MDEIT Jacobian for the same number of measurements, supporting the reconstruction analysis.

Implementation of a computational pipeline
Custom code in Matlab was integrated with the EIT solvers in EIDORS to solve the forward problem of MDEIT.A typical forward solution took approximately 10 min on a desktop computer.The Jacobian calculation was efficiently implemented using adjoint sensitivity analysis in COMSOL Multiphysics for which the Jacobian calculation took approximately five days to compute on a workstation computer.This is a substantial improvement on the forward method which was expected to require months of computational time and is the first time such a method has been implemented to the best of the authors' knowledge.The computational time could be further reduced to hours or minutes by writing and implementing a dedicated MDEIT Jacobian calculator, akin to EIDORS, which would sidestep the inefficiencies associated with COMSOL (Cheney et al 1990, COMSOL AB 2022).

Constant current versus constant voltage injection
A comparison of constant current and constant voltage injection showed no large difference in the SNR between them.For constant current injection, the MDEIT standing field will be caused by the current flowing in the wires and the head, and the change in magnetic field will be caused by redistribution of the current in the head since it can be assumed that the current does not redistribute within the wires.For constant voltage injection, the change in magnetic field will be caused by a combination of current redistribution and a change in the total current flowing through the wires and head.This means that the contribution to the signal from the wires can be considered to be negligible for constant current injection but not for constant voltage injection and can only be considered negligible for difference imaging.Incorporating the wires into the computational model or shielding the wires experimentally decreases the practicality of MDEIT, and modern EIT systems use constant current injection (Avery et al 2017).Therefore, constant current is preferable to constant voltage injection.

EIT versus MDEIT
1-axis MDEIT reconstructed images of higher quality than EIT in all cases.However, in some cases, the reconstruction quality was low for both techniques, (i.e.perturbation four and noise case one (figure 7(d))) and any comparison is of limited use.In general, for the lower SNR cases, EIT reconstructions tended to have a less clear boundary about the perturbation, with the amplitude of the reconstructed conductivity change decaying slowly as a function of the distance from the centre of mass of the perturbation whereas MDEIT reconstructions contained more discontinuities and were more fractured.When the SNR was larger as for perturbations one and two, MDEIT reconstructed images of superior quality to EIT regardless of whether one or three axes of the magnetic field were measured (figures 7(a), (b) and 9).1-axis MDEIT's ability to reconstruct superior images to EIT indicates that the enhanced quality cannot be attributed solely to the additional information obtained through tri-axial measurement, but rather suggests that an inherent advantage exists in MDEIT even when both methods had access to the same amount of information.EIT's greater rank deficiency can explain this difference, by showing that there was less independent information contained within the EIT Jacobian than the 1-axis MDEIT Jacobian.A physical interpretation of this is that the skull blurred the EIT signal more than the MDEIT signal, which matched expectations.The tank experiments presented here demonstrate that OPMs are feasible for use in MDEIT; however, this study cannot serve as a complete validation of fast neural MDEIT as it was conducted in two dimensions without a skull-like layer.A more robust validation would include image reconstruction in a three-dimensional headshaped tank with a skull-like layer.This was not practically achievable due to the magnetically shielded chamber available being too small for a realistic head-shaped tank and manipulation of an OPM in three dimensions.

Technical considerations
The EIT forward and Jacobian calculations were performed using EIDORS, whereas the MDEIT Jacobian was calculated in COMSOL Multiphysics, which necessitated an additional boundary and a simulated region of air around the human head.The forward model for MDEIT was calculated using EIDORS and custom code (with no such boundary conditions) (Polydorides andLionheart 2002, Polydorides et al 2022), meaning that there was less compatibility between the MDEIT solvers than for EIT.It was expected that this may decrease the quality of  the MDEIT reconstructions and favour EIT; however, this cannot be confirmed in this case until a fully integrated solver is created.
For the tank studies, the use of only one OPM in 25 positions instead of a full array of 25 OPMs was a technical necessity due to equipment limitations.It is expected that the movement of the OPM and the presence of a nearby stepper motor will have introduced more noise and uncertainty than could be achieved otherwise.For this reason, it can be expected that image quality would increase with a full array of OPMs and that real noise-based correction may perform better once these sources of uncertainty have been removed.

Future work
The implication of this is that there is rationale to proceed with performing simultaneous EIT and MDEIT measurements in a realistic tank and/or in vivo in a human with scalp electrodes and an array of magnetometers.Practically, the most suitable magnetometers for immediate use in MDEIT are SQUID magnetometers due to their large bandwidth of ∼1 MHz and high sensitivity of ∼1-10 fTHz −1/2 (Storm et al 2017, CTF 2021).The downside of using SQUIDs is that they are more expensive and less widely available or modifiable for future use than OPMs.For this reason, there is a rationale for working towards using OPMs instead of SQUIDs in the medium to long term.OPMs utilising helium-4 (MAG4Health 2021, Zahran et al 2022) have already been developed which have a larger bandwidth and dynamic range (∼2 kHz, ± 250 nT) than that of alkali vapour OPMs (∼100 Hz, ± 5 nT) which are popular for MEG (Shah and Wakai 2013, Hill et al 2020, Quspin 2022).This design is more compatible with the requirements of MDEIT but comes at the cost of a 4-fold decrease in sensitivity.Future work in the development of MDEIT should include an optimisation study for the design of OPMs for use with MDEIT.
It is possible that the image quality could be further increased by using 'magnetic injection, magnetic recording' rather than 'electric injection, magnetic recording' by inducing a current density in the brain with electromagnetic coils.This would then be MIT which has never been considered for neural imaging but is theoretically applicable (Watson andGriffiths 2001, Soleimani andLionheart 2006).MIT would eliminate the attenuation of the signal by the skull in both directions; however, there could be problems inducing the required current density at the optimal frequency range (∼1.5 kHz) and in the deepest structures of the brain.

Figure 1 .
Figure 1.The principle of MDEIT, showing the measurement of all three possible components of the magnetic field.Adapted from https://pixabay.com/illustrations/outline-head-profile-silhouette-5458467/.

Figure 2 .
Figure 2. A graphic representation of the FEM used in this work, showing the four perturbations considered.The perturbations are shown in black.All slices are taken through the centre of mass of the perturbation in each case.Adapted from Mason et al 2023.

Figure 3 .
Figure 3.A schematic of the setup used for the MDEIT tank experiments.

Figure 4 .
Figure 4.The convergence error ± mesh variability for the FEM convergence study.The error is expressed as a percentage of the standing field.

Figure 5 .
Figure 5.The SNR of the largest and mean top 10% of measurements for the magnetic field change due to four different perturbations for constant current and constant voltage injection.

Figure 6 .
Figure 6.The regularisation parameter versus the mean SNR of the largest 10% of changes in the magnetic field or voltage.

Figure 7 .
Figure 7.One example of image reconstruction using EIT, 1-axis MDEIT and 3-axis MDEIT for all noise cases and perturbations.Noise case 3 is empty for EIT reconstructions because noise cases 2 and 3 are identical for EIT.The slices are taken through the centre of mass of the true location of the perturbations.

Figure 8 .
Figure 8.The SNR of the single largest raw change and the mean of the largest 10% of changes as a function of the depth of the perturbation for all three noise cases.Adapted from Mason et al (2023).
4.5.MDEIT in a Saline tankMDEIT with an OPM successfully reconstructed clear images which corresponded to the true perturbation.TR with simulated NBC was the best algorithm, which matched previous expectations and is the algorithm used in fast neural EIT in vivo(Aristovich et al 2016, Faulkner et al 2018b).

Figure 9 .
Figure 9.The total reconstruction error (mean ± SE) for EIT, 1-axis MDEIT and 3-axis MDEIT for all noise cases and perturbations, showing the statistical significance between techniques for each perturbation (repeated measures ANOVA and multiple comparison test, N = 100).Adapted from Mason et al (2023).

Figure 10 .
Figure 10.The regularisation parameter for each of the 17 reconstructions with Tikhonov regularisation and NOSER.

Figure 11 .
Figure 11.One example of the unthresholded reconstructed perturbation for each reconstruction algorithm as well as the target perturbation.

Figure 12 .
Figure 12.The total reconstruction error for each reconstruction algorithm (mean ± SE) of the images after being thresholded at 50% the maximum negative change showing the statistical significance between techniques (repeated measures ANOVA and multiple comparison test, N = 17).The position error and area error composition of the total error are shown.
al 2018b).For N s samples, N m measurement channels and N e elements in the FEM, the input noise Î The two leading candidates for biomagnetic sensing of femtotesla scale fields are SQUID magnetometers and OPMs.Given the current state of technology, both technologies can meet the necessary bandwidth requirements, , Storm et al 2017lly having higher sensitivi, Brookes et al 2022, Zahran et al 2022scalp and allowing for a fully on-scalp system(Neuromag 2008, Storm et al 2017, CTF 2021, MAG4Health 2021, Brookes et al 2022, Zahran et al 2022, Cerca 2023).However, new OPMs are being developed with larger bandwidths (MAG4Health 2021), so it is possible that OPMs will be applicable to fast neural MDEIT in the future.1.5.7.Number of magnetometersThe number of magnetometers/measurement electrodes was chosen for four reasons.(1) Modern electroencephalography and magnetoencephalography data collection systems (such as those used in EIT) tend to have a number of channels that is a multiple of 16 (Brain Products 2016, Avery et al 2017) so 64 magnetometers/measurement electrodes satisfy this.

Table 1 .
Mason et al (2023)urbations considered in this work.The depth expresses the distance from the centre of mass of the perturbation to the surface of the brain.Reprinted fromMason et al (2023).

Table 3 .
Mason et al (2023)ses considered in this work for 1 and 232 measurement averages.Multiplicative noise is expressed as a percentage of the standing field.Adapted fromMason et al (2023).

Table 4 .
The number of rows, the rank and the ratio of the rank to the number of rows of the Jacobian matrix for EIT, 1-axis MDEIT and 3-axis MDEIT.