Skull Defects in Finite Element Head Models for Source Reconstruction from Magnetoencephalography Signals

Magnetoencephalography (MEG) signals are influenced by skull defects. However, there is a lack of evidence of this influence during source reconstruction. Our objectives are to characterize errors in source reconstruction from MEG signals due to ignoring skull defects and to assess the ability of an exact finite element head model to eliminate such errors. A detailed finite element model of the head of a rabbit used in a physical experiment was constructed from magnetic resonance and co-registered computer tomography imaging that differentiated nine tissue types. Sources of the MEG measurements above intact skull and above skull defects respectively were reconstructed using a finite element model with the intact skull and one incorporating the skull defects. The forward simulation of the MEG signals reproduced the experimentally observed characteristic magnitude and topography changes due to skull defects. Sources reconstructed from measured MEG signals above intact skull matched the known physical locations and orientations. Ignoring skull defects in the head model during reconstruction displaced sources under a skull defect away from that defect. Sources next to a defect were reoriented. When skull defects, with their physical conductivity, were incorporated in the head model, the location and orientation errors were mostly eliminated. The conductivity of the skull defect material non-uniformly modulated the influence on MEG signals. We propose concrete guidelines for taking into account conducting skull defects during MEG coil placement and modeling. Exact finite element head models can improve localization of brain function, specifically after surgery.

Magnetoencephalography (MEG) signals are influenced by skull defects. However, there is a lack of evidence of this influence during source reconstruction. Our objectives are to characterize errors in source reconstruction from MEG signals due to ignoring skull defects and to assess the ability of an exact finite element head model to eliminate such errors. A detailed finite element model of the head of a rabbit used in a physical experiment was constructed from magnetic resonance and co-registered computer tomography imaging that differentiated nine tissue types. Sources of the MEG measurements above intact skull and above skull defects respectively were reconstructed using a finite element model with the intact skull and one incorporating the skull defects. The forward simulation of the MEG signals reproduced the experimentally observed characteristic magnitude and topography changes due to skull defects. Sources reconstructed from measured MEG signals above intact skull matched the known physical locations and orientations. Ignoring skull defects in the head model during reconstruction displaced sources under a skull defect away from that defect. Sources next to a defect were reoriented. When skull defects, with their physical conductivity, were incorporated in the head model, the location and orientation errors were mostly eliminated. The conductivity of the skull defect material non-uniformly modulated the influence on MEG signals. We propose concrete guidelines for taking into account conducting skull defects during MEG coil placement and modeling. Exact finite element head models can improve localization of brain function, specifically after surgery.

INTRODUCTION
Localization of neuronal activity in the brain of patients is a common task in clinical neurophysiology. Specifically, during pre-surgical planning it is essential to locate physiological and pathophysiological brain activity as accurately as possible. Magnetoencephalography (MEG) is a non-invasive, functional recording modality that captures localizing information with high temporal resolution. Source reconstruction from MEG signals is a localization approach that reconstructs the distribution of neuronal currents inside the brain using a detailed volume conductor model of the patient's head.
Skull defects, such as post-surgical skull openings, are a challenge for functional evaluation of the brain. Initially, it was hypothesized that skull defects have a negligible influence on MEG signals and source reconstruction based on a small number of post-mortem phantom experiments, in vivo animal experiments, and simulation studies. Barth et al. (1986) used a physical coaxial dipole to simulate intracerebral currents in a formalin fixed human cranial specimen which was filled with conducting jelly. Okada et al. (1999b) recorded the MEG of a somatic evoked response in anesthetized piglets first over intact skull and then over the dura after a large skull section was removed (skull on vs. skull off). The main limitation of these experiments was that the skull defect was filled with nonconducting air instead of a conducting material mimicking the soft tissue in a healed skull defect. An early simulation of MEG above a 3 cm wide skull defect in a multi-sphere head model by Van den Broek et al. (1998) indicated that the MEG signals generated by sources placed next to a skull defect are not affected by the skull defect, and no source reconstruction errors could be observed for these source locations. However, sources under a skull defect remained an open question. Therefore, comprehensive evidence under realistic conditions that would allow one to generalize to post-surgical skull defects has been missing.
Recently, we reported on our in vivo animal experiment of the influence of conducting skull defects on MEG signals above and around two skull defects, using a well-defined current source under the middle and edge of one defect and next to it (Lau et al., 2014). The results demonstrated that skull defects can, in fact, substantially influence MEG signals and that the MEG signal magnitude deviates most if the source is under the middle of a skull defect. The change in MEG signal topography is dependent on the skull defect geometry and the relative position and orientation of the source. These measurements closed the gap in the experimental evidence. The question arises of whether volume conductor models of the head can reproduce these MEG signal changes accurately and whether source reconstruction from MEG in the presence of skull defects is possible.
The finite element (FE) method allows us to discretize the head volume into small volume elements and, therefore, to differentiate many tissue types and skull defects. Simulations by Vorwerk et al. (2014) showed that incorporating fine anatomical detail of a head with intact skull influences MEG signals and source reconstruction. Using this detailed modeling approach, Lew et al. (2013) investigated MEG in a FE model of a neonate head with and without differentiation of fontanels and sutures, which are natural skull defects. They reported relatively small differences in the MEG signals and source reconstruction for realistic source patches due to ignoring fontanels and sutures in the head model. The neonate skull is different from the adult one, because it is much thinner and of much higher conductivity. Consequently, the influence of non-neonate skull defects on MEG signals and source reconstruction remains to be modeled and evaluated. The MEG measurements of our controlled experiment (Lau et al., 2014) present a unique opportunity to validate a detailed FE model of skull defects with respect to corresponding in vivo MEG measurements and known source positions (Sander et al., 2010).
Therefore, the objectives of this study are (1) to assess the concordance of a detailed FE simulation of the controlled MEG skull defect experiment with the physical measurements of Lau et al. (2014), (2) to describe errors in source reconstruction from measured MEG signals caused by ignoring skull defects, (3) to assess the ability of an exact FE head model to overcome such source reconstruction errors, and (4) to identify critical modeling steps for skull defects.

Experimental Data
With ethics approval (Freistaat Thüringen, Germany, 02 034/08) we performed a 16-channel MEG of signals produced by a miniaturized artificial coaxial dipole implanted in a rabbit brain tangentially to the inner skull surface in vivo. The coaxial dipole had a contact at the tip that consisted of a 0.5 mm segment of exposed platinum wire of 0.25 mm diameter. The second ring shaped contact consisted of a 0.5 mm segment of the exposed platinum tube of 0.7 mm outer diameter. The distance between the contacts on the center axis was 1 mm, which resulted in a distance of approximately 1.1 mm between the contacts due to the conical shape. The dipole was connected to a constant-current source (20 Hz, 0.1 mA) with a sinusoidal waveform. Measured signals were band-pass filtered (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25) and approximately 300 consecutive trials (sinusoidal waves) were averaged. The resulting signal-to-noise ratio was above 40 dB, meaning that noise can be considered negligible in the analysis.
Following a reference recording with intact skull, two skull defects were introduced above the dipole and filled with conducting agar (1.0 S/m at 30 • C). The dipole was shifted inside the cortex under one of the skull defects (defect 1) in regular steps and measurements were taken at each step. The shifting was achieved by rotating a screw on the fixation device. One full rotation resulted in approximately 0.7 mm shift, called "full step, " and half of a rotation in a 0.35 mm shift, called "half step." We numbered the shift positions in 0.35 mm steps. Because the in vivo experimental setup had limited stability over time, we measured at full steps, i.e., every second shift position, except under one of the edges of defect 1 where the field-maps changed rapidly. To avoid tissue damage and bleeding that would compromise the results, the dipole implant was only moved in once to the skull defect site while the skull was intact.
Then, after introducing the skull defects, the dipole implant was shifted to the end of its range and stepwise back out under the defects. Consequently, intact skull recordings only exist for one side of defect 1, e.g., in Figures 2A, 5, 6. Further details of the experiment and measured signals are described in Lau et al. (2014).
The position of the artificial dipole was determined from postexperimental computed tomography (CT) with 0.4 mm isotropic voxels. The geometry of the skull defects was manually segmented from the CT. For each defect the outer edge, the inner edge and the middle of the cut surface was sampled finely to accurately represent its shape. We defined the shift eccentricity of a point on the shift line as 0 at the normal projection of the defect center on the shift line and +1 or −1 at the normal projections of the outermost edge points onto the shift line (see Figures 3, 5).
Before the experiment, a T2-weighted anatomical MRI data set (Siemens Magnetom 3T (Siemens, Erlangen, Germany), 224 slices, slice thickness 0.4 mm, 256 × 194 in-plane with resolution 0.423 mm, TR 2500 ms, TE 337 ms) was acquired of the anesthetized rabbit using a NORAS 8 Channel Multifunctional Coil (CPC) (NORAS MRI Products GmbH, Höchberg, Germany) that consisted of four coils directly above the rabbit's head and four coils below it. The signal intensity drop off with distance to the coils was corrected using nonparametric non-uniform intensity normalization (N3; Sled et al., 1998).
The co-registration was optimized to approximately 0.1-0.2 mm using a custom stereotactic frame, in which the rabbit head was fixated. Four spherical CT-visible markers were attached to the skull as landmarks to co-register the CT and head model. A well-defined configuration of coils, which was permanently mounted to the base of the stereotactic frame, was used to coregister the MEG sensors. The position of the rabbit head was continuously monitored.

Segmentation and Meshing
The tissues of the head were segmented semi-automatically using Seg3D 1 from the co-registered CT and MRI volume data sets. For computational efficiency, the ears, the neck, and the lower part of the head, consisting primarily of the jaw, were removed from the body compartment. The skull was segmented by thresholding the CT and manually correcting for the artifact of the dipole implant and discontinuities due to the limited voxel resolution. Natural skull openings toward the spine and at the optic nerve exit were modeled open. The burr hole at the posterior lateral region of the skull, through which the coaxial dipole was inserted, was closed in the segmentation, because it was sealed with plastic and nonconducting glue during the experiment. It was ensured that the cancellous bone was enclosed by compact bone.
The white matter, CSF, and major blood vessels were segmented using thresholding of the T2 volume data set and manual completion using an anatomical atlas of the rabbit brain. The CSF layer around the brain was much thinner than one voxel side length in the pre-experimental MR image and, during the experiment, the brain was pressing against the intact, almost transparent dura at the interior boundary of the skull defect. Consequently, no CSF layer around the brain was differentiated.
The skull defect geometries were manually segmented from the CTs as accurately as the voxel resolution allowed. The liquid layer was derived by dilating the skull by 1 voxel and subtracting the body (including the skull). The stem of the coaxial implant was not differentiated. The miniaturized coaxial design eliminated most interference.
The segmentation was converted into a hexahedral mesh from the voxel coordinates using the software Vgrid version 1.3.1 2 as shown in Figure 1. Using Vgrid, node-shifting (Camacho et al., 1997;Wolters et al., 2007) with a smoothing factor of 0.49 was applied to smooth the surfaces between compartments.

Conductivities
Equivalent homogeneous conductivity values, summarized in Table 1, were assigned to each segmented tissue type as well as the agar in the skull defects and a layer of conducting liquid on the skull surface mimicking a thin layer of skin. These conductivity values were used for both forward simulations and source reconstructions. The tissue conductivity values were taken from the literature (Geddes and Baker, 1967). The scalp/body conductivity was chosen as a mixture of fat, skin, and muscle, weighting fat strongest due to the subcutaneous fat layer (Gabriel et al., 2009). For the inner and outer compact bone layers of the skull a conductivity of 0.004 S/m (Tang et al., 2008) was used. Because the cancellous bone consists of bone, fat, cells, and blood, an approximate value of 0.046 S/m (Haueisen et al., 1997;Akhtari et al., 2002Akhtari et al., , 2006 was used. Gray matter conductivity of 0.23 S/m was used (Crile et al., 1922;Akhtari et al., 2006). For white matter, an equivalent isotropic conductivity of 0.31 S/m was chosen, which lies between the lower conductivity perpendicular  to the fibers (0.125 S/m measured by Nicholson, 1965) and the higher conductivity parallel to the fibers (1.18 S/m measured by Nicholson, 1965). The conductivity of the vitreous humor of the eye was measured to be 1.55 S/m at 25 • C (Pauly and Schwan, 1964). Since the eyes of a rabbit are positioned mostly outside the skull, a temperature of 30 • C was assumed and a conductivity of 1.55 S/m + 5 * 0.031 S/m = 1.7 S/m was used for the vitreous humor and the aqueous humor. The conductivity of the lens was estimated with 0.35 S/m based on the measurements of Pauly and Schwan (1964) and Lindenblatt and Silny (2001).

MEG Coil Geometries and Positioning
The geometric setup of the in vivo rabbit experiment (Lau et al., 2014) was reproduced in the simulation setup as far as possible. The exact asymmetric MicroSQUID gradiometer array geometry (Nowak et al., 1999) was used, except for the base length of the gradiometers that had to be modeled uniformly with 30 mm instead of 28.5-31.5 mm to meet the requirements of the simulation software. The field-map magnitude and topography deviation between simulated nonuniform and uniform base length geometries in the experimental setup was small (MAG rel <0.001, RDM * <0.01).

Forward and Inverse Solution Methods
Various approaches have been developed to solve the MEG forward problem and model the source singularity. Based on recent comparison studies with regard to numerical and modeling accuracy as well as computational speed (Wolters et al., 2007;Lew et al., 2009;Vorwerk et al., 2014), we used the Venant direct FEM approach for modeling the dipole source (Buchner et al., 1997;Wolters et al., 2007) and trilinear basis functions in an isoparametric FEM approach in the geometryadapted hexahedral volume conductor model (Wolters et al., 2007). This approach has a high computational efficiency when used in combination with the FE transfer matrix approach (Wolters et al., 2004) and an algebraic multigrid preconditioned conjugate gradient solver (Lew et al., 2009). We performed our computations using the SimBio software toolbox. 3 In forward simulations, source strength of 130 µAmm was used, given that a constant current of 100 µA was injected, that the shortest distance between the electrodes was approximately 1.1 mm, and that the current spreads across the dipole contact surfaces to some degree. Given the a priori knowledge of a single dipolar source, an unconstrained moving dipole fit was chosen to reconstruct the source properties from the in-vivo MEG measurements. To speed up convergence, the source position derived from the CT was used as the initial guess. Downhill simplex optimization (Nelder and Mead, 1965) was used (reflection factor −1.0, expansion factor 2.0, contraction factor 0.5, max. 200 iterations). The initial simplex size of 1 cm was chosen to match the dimensions of the rabbit brain, which is approximately 3 cm wide from left to right.
A limitation of the MicroSQUID gradiometer array was that it is mono-directional. This means that source components that are oriented normal to the MEG coil plane are not well detected (Haueisen et al., 2012). Consequently, after source reconstruction, the source orientation component that is normal to the MEG coil plane is not well determined by the measurement data of this system. The experiment was designed so that the implanted source was oriented almost parallel to the MEG coil plane (angle approximately 3 • ). Therefore, orientation components normal to the MEG coil plane were minimal and could be excluded from the analysis by setting this orientation component of the equivalent source to zero.

Configurations for Source Reconstruction
During source reconstruction, three types of configurations were used (see Figures 5, 6):

Reference (i-I): The intact skull measurements (denoted i)
were reconstructed using the intact skull FE model (denoted I). This configuration is color-coded in green. 2. Ignoring skull defects (d-I): The measurements in the presence of skull defects (denoted d) were reconstructed using the intact skull FE model. This configuration is color-coded in red. 3. Incorporating skull defects (d-D): The measurements in the presence of skull defects were reconstructed using the FE model incorporating the skull defects. This configuration is color-coded in blue.

Field-Map Measures
To quantify the topographical deviation caused by defect 1 and by the combination of both skull defects, we determined the relative difference measure RDM * (Meijs et al., 1989) where i and j are the channel indices, v i,r is the value of the reference signal with intact skull, and v i,s is the value of the signal measured with either one or two defects. The RDM * is a goodness of fit measure that indicates the difference of two multichannel data vectors, such as two MEG field maps. Both sample vector values are divided by the respective vector norm to exclude scaling differences, which relate to the source strength rather than location and orientation. The RDM * ranges from 0 to 2, where 0 means that the multichannel data vectors are identical and 2 means that one is the inversion of the other. To quantify the magnitude deviation caused by these two conditions, we determined the relative magnitude difference MAG rel (Güllmar et al., 2006): For the comparison of simulated intact skull field maps with simulated skull defect field maps in the shift series, the reference map was the respective intact skull field map with the source at the same position as in the skull defect field map being evaluated.
In the in vivo rabbit experiment, intact skull field maps could only be acquired for few of the shift positions. Therefore, the reference map for the comparison of measured field maps was the intact skull field map obtained with the dipole positioned close to the center of the subsequently introduced defect 1.

Equivalent Source Characteristics
The reconstructed positions of the sources were characterized with the Euclidean distance to the physical position during the experiment. Taking into account the radius of the outer contact of 0.4 mm, the expected ideal value range for the Euclidean distance was assumed to be 0 to 0.4 mm. The spacing of equivalent sources along the shift line was characterized with the distance to the equivalent source one full step along the shift line. The expected ideal value is the physical distance of sources of 0.7 mm. The angle between the equivalent source orientation and the physical source orientation was evaluated. Angle variations of up to 16 degrees were within the expected ideal range, because the conic shape of the coaxial dipole surface had an angle of 16 • between the center line and the line connecting the inner and outer platinum contact. The strengths of the equivalent sources were evaluated relative to an expected ideal range of 100-130 µAmm, given that a constant volume current of 100 µA was flowing between the platinum contacts that were approximately 1.1 mm apart at their closest points. The explained variance of the reconstruction was reported. Two groups of physical source positions were selected for comparison (see, e.g., Figure 5): (1) source positions under defect 1 (shift eccentricity −0.5 to +0.5, shift position 15-22) were marked with a square, (2) source positions next to defect 1 and 2 (shift position 6-10) were marked with a sphere.

Sensitivity Analysis
To investigate the sensitivity of the source reconstruction to tissue conductivities, the conductivities of compact bone, cancellous bone, gray matter, and white matter were varied by ±10, ±20, and ±50% ( Table 2).

Forward Simulations vs. In vivo Measurements
The MEG signals were simulated forward for a series of source positions along a linear path under defect 1, as shown in Figure 2D. These source positions matched the physical positions of the implanted artificial coaxial dipole during the in vivo measurements ( Figure 2C). The simulations were repeated with the intact skull model ( Figure 2B) to evaluate their concordance with the available intact skull recordings (Figure 2A).
The measured and simulated MEG signals above intact skull are shown side by side in Figures 2A,B for qualitative comparison. The gradual change in the measured flux density map topography as the source was shifted from one side of defect 1 to the other was reflected clearly in the simulated flux density maps ( Figure 2C vs. Figure 2D). The absolute magnitude of the simulated flux density maps with constant source strength of 130 µAmm was within approximately ±30% of that of the recorded flux density maps (Figure 2).
Figures 3A,C shows the magnitude and topography differences between intact skull MEG signals and MEG signals above skull defects in measurement and simulation. In the measurements of the rabbit, the MEG signal magnitude dropped by approximately MAG rel = 0.2 (20%) when the source was central under defect 1 (Figure 3A, blue square markers). The simulated MEG signals ( Figure 3A, blue 1.0 S/m curve) showed the same magnitude reduction of approximately MAG rel = 0.2 for the source central under defect 1.
The simulated topographic difference between the intact skull flux density maps and skull defect flux density maps (Figure 3C, 1.0 S/m) was approximately RDM * = 0.07 for sources under defect 1. In the measurements of the rabbit, this difference was a bit larger (RDM * approximately 0.16), which could be due to experimental variability. The slight discrepancy between measurement and simulation of RDM * values for MEG at the lowest and highest shift positions is partly due to the fact that for the measured data one intact skull measurement (Figures 3A,C, diamond marker) close to the center of defect 1 was chosen as the reference from the limited set of available ones. In the simulation, intact skull MEG signals could be computed and used as the respective reference for all source positions. Therefore, the forward  . The dipolar source is indicated with a black bar with two spheres marking the poles. Skull defects are marked by closed black lines indicating the inner, middle, and outer boundaries of the defects (see Figure 5 for a three-dimensional view). MEG coil positions are marked with gray dots. The minimum and maximum value (in fT) and the iso flux density line step are displayed above each map.
simulations reflect best the influence of skull defects on the magnitude and topography at the lowest and highest shift positions.
To evaluate the influence of the skull defect conductivity on the MEG signal change, the forward simulations were repeated with skull defect conductivities of 10 −6 S/m (approximately  (Figure 3, color-coded). For skull defect conductivities below that of intact three-layer skull, i.e., 10 −6 S/m and 0.004 S/m, the MEG signal magnitude was marginally increased (MAG rel = 0.02) for a source under a skull defect. For defect conductivities above that of intact three-layer skull, the MEG signal magnitude was decreased. In Figures 3B,D, the simulated magnitude and topography differences due to the presence of the skull defects are plotted over the range of defect conductivity values. The MEG signal magnitude differences increased as the defect conductivity was increased and this effect was greater for small increases of the defect conductivity above that of skull ( Figure 3B). The MEG signal magnitude change in this setup reached 50% of the largest observed magnitude change per source position in a range of approximately 0.2-0.3 S/m defect conductivity, which emphasized this non-uniformity. The flux density map topography change showed a similar non-uniform dependency and reached 50% of the largest observed topography change per source position also in a range of approximately 0.2-0.3 S/m ( Figure 3D).
To investigate the nature of the topography change in the forward simulations, the intact skull flux density maps were subtracted from the skull defect flux density maps at different source positions, which are shown in

Intact skull measurements with intact skull model (i-I)
The intact skull equivalent dipoles (

Skull defect measurements with intact skull model (d-I)
The equivalent dipoles of the skull defect MEG measurements reconstructed using the intact skull FE model (Figure 5, red) were arranged in sequence along the shift line but not equidistant. The sources next to the defects (round markers) had an offset from the shift line of approximately 0.51 mm [0.43 mm, 0.53 mm], which is comparable to the i-I result. Sources under defect 1 (square markers), however, were displaced radially away from defect 1, partly to the side of defect 1 (Figures 5A,C) and partly deeper ( Figure 5C)

Skull defect measurements with skull defect model (d-D)
When the exact FE head model incorporating the skull defects was used to reconstruct the skull defect MEG signals (Figure 5 To investigate the global offset of equivalent sources from the shift line, all equivalent sources were plotted relative to the physical dipole at the time of recording (Figure 5B). The d-I reconstructions of sources under defect 1 (red squares) were displaced from the implant surface, while the remaining sources were arranged close to the conic implant surface. Defect Conductivity Variation Figure 6 shows the influence of the defect conductivity in the model on the equivalent sources. For defect conductivity 0.004 S/m, the equivalent sources reconstructed using the skull defect model (d-D reconstruction, blue) were positioned and oriented most similarly to the intact three-layer skull model (d-I reconstruction, red) (Figure 6). For defect conductivity 10 −6 S/m (quasi-non-conducting), the displacement of equivalent sources was stronger (Figure 6, column 1). As the defect conductivity increased toward the experimental value of 1.0 S/m, the equivalent sources under defect 1 were positioned increasingly on a line (Figure 6). The distance to the physical positions (Figure 6) of equivalent sources diminished as the defect conductivity of the model approached that of the experiment, with the greatest change at the initial conductivity steps (0.4 S/m and below). The orientations of equivalent sources next to the defects approached the shift direction as the defect conductivity increased, with the greatest change for small increases of the defect conductivity. Equivalent sources under defect 1 were only marginally reoriented. The strength of the equivalent sources non-uniformly increased as the defect conductivity increased.

Concordance of Measurement and Simulation
The forward simulation of the flux density maps reflects all topographic and magnitude features that were observed in the in vivo animal experiment. In particular, the gradually changing flux density map difference between the intact skull and the skull defect condition observed in the in vivo experiment (Figure 8 in Lau et al. (2014)) was also observed in the FE simulation (Figure 4). This agreement confirms that the observed changes pertain to the skull defects and are not due to experimental or modeling limitations. The conducting medium occupying the skull defect volumes represents a conducting bridge into and through the skull. This causes a displacement of the volume currents into and through the defect.
The magnitude difference between the measured and simulated flux density maps in the intact skull condition is likely due to experimental variability of the tissues around the implanted source. The intact skull recordings at the start of the experiment required approximately 15% higher voltage than the skull defect recordings several hours later to produce the same current. This means that the net resistance around the source was higher, i.e., the net conductivity was lower. However, the current remained constant.
In agreement with the in vivo experiment, the simulated flux density maps above conducting skull defects were approximately 20% weaker than the intact skull equivalent for sources that were central under defect 1. The magnitude reduction was diminished toward the boundaries of defect 1 in both measurement and simulation ( Figure 3A). This confirms that a detailed FE model with conducting skull defects can quantitatively represent the influence of skull defects. A FE simulation of a neonate head (Lew et al., 2013) differentiating a single-layer non-ossified skull (conductivity 0.04 S/m) from fontanels and sutures (conductivity 0.3 S/m) also produced a change of MEG signal magnitude of up to absolute MAG rel = 7.9% due to the presence of the fontanels and sutures. Considering that the neonate skull was very thin and had comparatively high homogeneous conductivity, and that the absolute MAG rel was calculated over the whole head, including large field-map areas distant from defect 1, these absolute values are in accordance with our results.
The conductivity of the skull defects has a non-uniform effect on the flux density map change due to skull defects ( Figure 3B). Small increases of the defect conductivity above that of intact three-layer skull have the strongest effect. The presence of a higher conducting path through the weakly conducting skull is the dominant reason for the signal change. This was described as the shunting effect (Brody, 1956;Rush and Driscoll, 1968;Haueisen et al., 1997). The MEG signal magnitude and topography change in this setup reached 50% of largest observed magnitude and topography change per source position at defect conductivities of 0.2-0.3 S/m, which is in the range of soft tissues.
Non-conducting skull defects, such as air-filled ones during surgery, have comparatively little influence on the flux density map. Our simulation of skull defects with conductivity 10 −6 S/m caused a flux density map magnitude increase of less than 2%. Similarly to intact skull, the non-conducting defects represent a barrier to the volume current. When the defect conductivity approaches zero, the small amount of volume current that could have passed through the skull is also displaced into the brain volume. This is consistent with earlier observations involving non-conducting skull defects (Barth et al., 1986;Okada et al., 1999a,b).
An early simulation of human MEG using a spherical finite element head model (Van den Broek et al., 1998) investigated the flux density map topography above a skull defect of 3 cm diameter generated by tangential sources next to the skull defect at different depths and found no topography change. Their source positions are approximately equivalent to shift eccentricity ±1.3 (next to skull defect) in our study, where the influence of the skull defect is small in our study as well.

Intact Skull FE Model vs. Skull Defect FE Model
The reconstruction from the intact skull flux density maps (i-I configuration) serves as a reference condition. The full step distance of equivalent sources of 0.7 mm matches the physical distance of 0.69 mm very well. This indicates internal consistency of the model. The equivalent sources are arranged at the boundary of the conic surface of the implant tip ( Figure 5B) with an effective distance to the shift line that is comparable to the outer radius of approximately 0.4 mm of the larger platinum contact. This indicates that the volume current was not completely symmetric around the shift line due to differences in the conductivity of the tissues surrounding the implant. The observed small angle between shift line and equivalent source orientations of approximately 6 • can then be understood as a consequence of the conic shape of the physical source (angle of 16 • to the shift direction).
Ignoring conducting skull defects in a head model (d-I configuration) can cause errors in position, orientation and strength of equivalent sources reconstructed from MEG signals. Sources under a skull defect are mostly displaced away from the defect; partly toward the boundary of the defect and partly deeper (Figures 5A,C). This matches the reduced MEG signal magnitude above the defect (Lau et al., 2014). When the skull defect is ignored, the source should be estimated farther away from the defect. In this case, the equivalent sources were displaced under or next to the boundary of the skull defect ( Figure 5A). This displacement might scale to larger skull defects. The larger a skull defect, the more dominant is its influence on the flux density map and consequently on the equivalent sources. Sources next to a skull defect are mostly reoriented; in this experiment by approximately 35 • . This is explained by the reduced MEG signal magnitude above the defect causing a net reorientation in the flux density map. Such orientation errors should be observable in sources that are proximal to skull defects in humans. The strength of equivalent sources next to the defects was underestimated by approximately half due to ignoring the skull defects. Van den Broek et al. (1998), in their simulation study, used a spherical finite element head model with an intact skull layer to reconstruct sources next to a skull defect (equivalent to shift eccentricity ±1.3), which was present during forward simulation of the MEG signals. In agreement with our study, they found no location errors. The fact that they did not observe orientation errors is likely due to the complete symmetry in their model. The volume conductor was spherical, the defect had a circular shape, the source orientation was in plane with the center of the defect and the source was exactly tangential to the spherical skull surface.
A FE simulation study (Lew et al., 2013) of a neonate head with fontanels and sutures also produced position, orientation and magnitude errors in source reconstruction from MEG signals due to ignoring the fontanels and sutures. However, because of the much smaller conductivity difference between the thin non-ossified skull and the fontanels and sutures, their errors were smaller relative to the skull defect sizes.
Modeling skull defects in a realistic FE head model can compensate the errors in position, orientation and strength of equivalent sources. In our study, the displacement of equivalent sources away from defect 1 was corrected to reflect the linear sequential arrangement of physical sources. A small spatial dispersion of sources under defect 1 remained, which indicates the mesh resolution of the model may be at the limit that is necessary to exactly represent the fine defect compartment boundary. The mesh size was derived from highest available CT imaging and MRI resolution of 0.4 mm 3 . The benefit of a smaller element size derived from even higher resolution volume imaging should be investigated further. The orientation error of sources that are physically next to the defects was fully compensated by the incorporation of the skull defects in the model (Figure 5A, blue dots). The strength of these sources was restored from approximately half (41 µAmm) of the i-I reference strength (78 µAmm) to a value very close to it (80 µAmm).

Role of Defect Conductivity
The modeled defect conductivity non-uniformly modulates the influence of a skull defect on the position, orientation and strength of sources reconstructed from MEG signals. Small increases of defect conductivity above that of compact bone in the model most strongly reduce the position, orientation and strength error. This corresponds to the finding that the forward simulated flux density map changes most strongly for small increases of the defect conductivity above that of compact bone (Figure 3). Already for small defect conductivity values, e.g., below 0.2 S/m, the forward simulation explains most of the measured flux density map. In humans, non-acute skull defects would be occupied by tissue with conductivity around or below 0.2 S/m. Consequently, the demonstrated source reconstruction errors may be observed in humans. In contrast to the orientation and strength error, the position error of sources under a defect only gradually reduces as the defect conductivity in the model approaches the physical one (Figure 6). For accurate localization of sources under skull defects, the modeled defect conductivity should be chosen to match the physical defect conductivity.

Sensitivity to Tissue Conductivities
The variation of compact bone, cancellous bone, gray matter, and white matter tissue conductivities (not shown due to volume) confirmed that the values derived from measurement literature (Section Conductivities) are appropriate and that the source reconstruction results from MEG signals are stable in an environment of at least ±10% around the individual conductivity values. Larger changes of the conductivity of an individual tissue type gradually displaced the overall set of equivalent sources along the superior-inferior axis by small amounts. Gençer and Acar (2004), using their simulated voxel-wise maps of the sensitivity of the MEG signals to tissue conductivity changes, found that MEG signals, and with them their source reconstruction, are more sensitive to changes of the conductivity of the skull in the vicinity of the source. This is in accordance with our results.
The results indicate that strong changes to tissue conductivities of ±20% or more can influence source reconstruction from MEG signals and that tissue conductivities should be chosen as accurate as possible, which is in agreement with a sensitivity study of Van Uitert et al. (2004) using a human head shape. A small variation of the compact bone conductivity had comparatively little influence on source positions, orientations and strengths, which was also found by Stenroos et al. (2014). This further supports the observation that the origin of the source reconstruction errors is the presence of a highly conducting path through a weakly conducting compact bone layer, rather than the particular conductivity values. The influence of cancellous bone conductivity on position is stronger for a model with skull defect than a model without skull defects. This indicates that the skull defect allows the volume current to pass through the cancellous bone layer. Consequently, it is important to differentiate cancellous and compact bone layers in head models that involve post-surgical skull defects.

Modeling Skull Defects in Humans
The source reconstruction errors due to ignoring skull defects in adult human MEG investigations are expected to be substantial. Human post-surgical skull defects are usually larger than in this experiment, in absolute size as well as relative to the skull surface. Increasing the skull defect size is known to cause a more substantial displacement of volume currents, a wider spread of electroencephalography (EEG) signal changes (Van Burik and Peters, 2000;Oostenveld and Oostendorp, 2002) and larger EEG signal source reconstruction errors (Lanfer et al., 2012). The same principle applies to MEG signals. The increased displacement of volume currents causes wider spread MEG signal changes, which result in larger source reconstruction errors. The neonate human skull should be differentiated from the adult one, because it is much thinner and of higher conductivity. The fontanels and open sutures present a weaker conductivity contrast to the non-ossified skull and consequently their influence should be weaker than in adults. This is supported by an EEG source reconstruction study in five neonates (Roche-Labarbe et al., 2008) as well as a MEG-EEG finite element simulation using a neonate head (Lew et al., 2013).
Sources close to a skull defect need to be differentiated from distant sources. Ignoring skull defects during source reconstruction from MEG signals displaces sources that are physically under a defect away from that defect. Sources physically next to a skull defect are reoriented and reduced in source strength. While in humans, the types of source reconstruction errors are the same, the much larger skull defect volume and large head surface portion of a skull defect is expected to cause larger position, orientation and strength errors than in this study. These errors occur at defect conductivities in the range of human soft tissue. Consequently, natural skull openings, such as the orbital fissures, should also be included in head models. The conductivity of post-surgical skull defect tissue over time should be investigated further, because it has a non-uniform influence on MEG signals. Further, the influence of model errors and skull defects, respectively, on the confidence region around the optimal source location and orientation should be investigated in human head geometry with whole-head sensor coverage.
The direct comparison of MEG measurements of a controlled source at known locations under a skull defect with an exact finite element forward simulation of the same setup validates the finite element approach to volume conductor head modeling in the presence of skull defects qualitatively and quantitatively. An alternative approach to head modeling is the boundary element (BE) method, which is based on a discretization of tissue compartment boundaries. The standard BE approach had been to assume nested compartments, such as scalp, skull, and brain (Kybic et al., 2006). This partly allowed for modeling skull defects in EEG (Bénar and Gotman, 2002;Oostenveld and Oostendorp, 2002) either by joining the scalp, defect, and brain compartment (Van Burik and Peters, 2000), by thinning the skull compartment at a skull defect (Roche-Labarbe et al., 2008), or by introducing an extra defect compartment inside the skull layer. Recent developments of the BE method relaxed the limitation of nested compartments for EEG (Kybic et al., 2006;Stenroos and Sarvas, 2012) and MEG (Stenroos and Sarvas, 2012). The application of BE approaches to skull defects in MEG should be investigated further and validated.
The spatial accuracy of source reconstruction from human MEG signals is determined by a combination of factors, including signal to noise ratio, spatial sampling density, co-registration accuracy, sensor coverage, and head model simplifications such as ignoring skull defects. In our controlled experiment, we isolated the influence of skull defects by minimizing the impact of all other factors. We showed that realistically incorporating skull defects in head models enables us to reconstruct sources from MEG signals recorded in the presence of skull defects. In a human application setting, this approach reduces the compound error by eliminating one error source. MEG studies should account for conducting skull defects during pick-up coil placement and modeling in the following way: 1. The MEG signal should be sampled densely in the proximity of a skull defect to capture the higher spatial frequencies. 2. The co-registration accuracy of MEG coil positions relative to the head should be maximized. 3. The skull defect boundary should be segmented as exactly as possible, if available from a CT, because CT is not geometrically distorted and has a good bone-tissue contrast. 4. The types of tissue occupying the defect volume should be determined and differentiated carefully and matching conductivity values should be used. 5. The compact and cancellous bone layers adjacent to skull defects should be differentiated and modeled as exactly as possible, because of the leakage of current into the cancellous bone.

CONCLUSION
Our study provides qualitative and quantitative evidence based on experimental measurements that ignoring skull defects in volume conductor head models can influence source reconstruction from MEG in terms of location, orientation, and strength of the equivalent sources. The skull defect conductivity has a non-uniformly modulating influence on the changes to signal and equivalent source. A detailed realistic finite element model incorporating skull defects can reproduce the influence of skull defects on the MEG signals. Realistic FE head modeling combined with appropriate instrumentation enables source reconstruction from MEG in the presence of skull defects.
We conclude that skull defects should be incorporated in volume conductor models of the head for source reconstruction from MEG signals. This is particularly important in postsurgical and post-traumatic cases in which the brain area under investigation is typically under or proximal to the skull defect. The conductivity of skull defect tissues should be investigated further. Our next step is to investigate the utility of exact FE head models in source reconstruction from EEG in the presence of skull defects using the same controlled-source experiment. This modeling approach can improve the diagnostic localization and characterization of brain activity, such as epileptic discharges, in post-surgical patients.

AUTHOR CONTRIBUTIONS
SL and JH conceived the study. SL designed and developed the experimental setup and experimental methods. SL performed the data acquisition. DG facilitated the MR imaging of the animals. LF facilitated the anesthesia and surgical work on animals. SL processed the measured data and volume data sets. SL constructed the FE model and performed the simulations and source reconstructions. SL performed the analysis, produced the images, interpreted the results, and wrote the paper. All authors took part in the scientific discussion at multiple stages of the study and provided feedback from the experimental (JH, LF), modeling (JH, CW, DG), theoretical (JH, CW, DBG), and clinical (MC, DBG) perspective. All authors reviewed the manuscript and approved it for publication.