Balanced, bi-planar magnetic field and field gradient coils for field compensation in wearable magnetoencephalography

To allow wearable magnetoencephalography (MEG) recordings to be made on unconstrained subjects the spatially inhomogeneous remnant magnetic field inside the magnetically shielded room (MSR) must be nulled. Previously, a large bi-planar coil system which produces uniform fields and field gradients was used for this purpose. Its construction presented a significant challenge, six distinct coils were wound on two 1.6 × 1.6 m2 planes. Here, we exploit shared coil symmetries to produce coils simultaneously optimised to generate homogenous fields and gradients. We show nulling performance comparable to that of a six-coil system is achieved with this three-coil system, decreasing the strongest field component Bx by a factor of 53, and the strongest gradient dBx/dz by a factor of 7. To allow the coils to be used in environments with temporally-varying magnetic interference a dynamic nulling system was developed with a shielding factor of 40 dB at 0.01 Hz. Reducing the number of coils required and incorporating dynamic nulling should allow for greater take-up of this technology. Interactions of the coils with the high-permeability walls of the MSR were investigated using a method of images approach. Simulations show a degrading of field uniformity which was broadly consistent with measured values. These effects should be incorporated into future designs.

Magnetoencephalography (MEG) non-invasively measures magnetic fields outside the head generated by neuronal currents in the brain, providing information on brain function with high spatial and temporal resolution 1,2 . The MEG signal is commonly on the level of tens of femtotesla (fT), significantly weaker than the magnetic interference from many other sources 2 . To measure these minute fields, state-of-the-art MEG systems use superconducting quantum interference devices (SQUIDs) and are housed inside a magnetically shielded room (MSR). The SQUIDs generally require cooling to a temperature of 4 K using liquid helium. The need for cryogenic cooling results in a one-size-fits-all sensor arrangement where the subject must remain with their head fixed in position inside the MEG helmet during a scan. This limits the range of subject groups that can be investigated, since, for example, a child's head will be far from the sensors, resulting in decreased sensitivity to brain activity, and although baby-focused MEG devices do exist, lifetime-compliance cannot be achieved with a single system 3 . The neuroscientific questions that can be addressed are also restricted, as subjects have to remain relatively still during measurements. Implementation of tasks requiring significant subject movement, such as social interaction or spatial navigation is therefore challenging.
As a result of the limitations of SQUID-based measurements of neuromagnetic fields, MEG systems built using optically pumped magnetometers (OPMs) are garnering considerable interest [4][5][6][7][8][9] . OPMs measure changes in magnetic field by exploiting the optical properties of a heated vapour cell containing spin-polarised alkali atoms 10 . As a result of a considerable effort in miniaturisation 6 , it is now possible to produce OPM sensors which can be flexibly placed on the scalp of a subject. Simulation studies suggest that a four-fold increase in signal strength and improved accuracy in source reconstruction can be achieved through the resulting reduction in source-to-sensor distance relative to typical SQUID-based MEG systems 11,12 . With such a wearable setup the subject is also notionally free to move during an experiment, provided the OPMs remain in a fixed position with respect to the head.
QuSpin OPMs are commercially available magnetic field sensors, which operate in the spin exchange relaxation free regime with a sensitivity of <15 fT/√Hz 13 . In this manuscript we make use of two sensor types: Generation 1 and 2 QuSpin sensors. Generation 1 sensors have a sensitive volume of gas located at ~6 mm from the sensor surface with external dimensions 13 × 19 × 110 mm 3 . Generation 2 sensors have an approximately 6.5 mm distance from the sensitive volume to outer surface, but with a more compact casing of 12 × 17 × 24 mm 3 size (https://quspin.com/products-qzfm/). With these sensors, high sensitivity operation can only be achieved for fields measured within a narrow operational range of up to ±5 nT 13 . The passive shielding employed in a typical 2-layer mu-metal MSR, reduces the remnant Earth's field to a few 10's of nT. Further reduction of the field in the vapour cell can be achieved by using "on-sensor" coils wrapped around the vapour cell. At the start of an experiment these coils produce a compensation field, locked and optimised for a fixed sensor position and orientation 14 . However, with this approach translation or rotation of the OPM in the spatially varying remnant field inside a MSR can easily produce field changes which exceed the sensor's operational range, thus rendering the sensor's output unusable until it returns to its original position.
Our previous work 5,7,8 showed that the remnant field in a MSR can also be nulled using a set of six, bi-planar coils mounted on large planes of dimensions 1.6 × 1.6 m 2 , separated by 1.5 m, which compensate B x , B y and B z as well as the three dominant first-order spatial field gradients G x = dB x /dz, G y = dB y /dz and G z = dB z /dz. These coils were used in a MSR (Vacuumschmelze, Hanau, Germany) comprised of two layers of mu-metal and one layer of aluminium to reduce the dominant static field component, B x , by a factor of 46 from 21.8 ± 0.2 nT to 0.47 ± 0.08 nT, and the dominant gradient, G x , by a factor of 13 from 7.4 nT/m to 0.55 nT/m 7 . The integration of these coils with a scalp-mounted OPM array allowed high quality MEG data to be recorded whilst the subject made large (±10 cm and ±35°) movements 7 . With this set-up, MEG recordings have been made while unconstrained subjects perform naturalistic experiments, such as a 'bat and ball' game 5 , as well as when performing a language lateralisation task 8 .
The bi-planar coils thus form a crucial element of an OPM-MEG scanner that allows measurements to be made on a moving subject, but their construction is time-consuming due to the large size of the coils and the significant number of coil layers that are required. The distribution of the weight of the current carrying wires and structural rigidity of the coil-set also have to be taken into account during coil construction.
Here, we show that by exploiting the natural symmetry shared by pairs of coils (B x and G x , B y and G y and B z and G z ) the number of coils that need to be constructed to produce a simple field-nulling system can be reduced from six to three, as illustrated in Fig. 1. We describe the harmonic-minimisation approach that allows optimisation of planar current distributions to produce both homogenous fields and spatial field gradients. The performance of coils generated using the joint design method is compared with that of coils designed to produce only a uniform field or field gradient. We also detail a proportional integral derivative (PID) control method, which Each coil pair is simultaneously designed to produce a homogeneous field and a homogeneous field gradient over a region at the centre of the coils. The three coil layers on each plane are spaced by ~1 mm in the constructed system. nulls the static field and field gradients automatically by independently optimising the currents in each plane based on reference measurements. The static field compensation achieved using this method is quantified by comparing field maps measured in a single plane in the MSR before and after field compensation.
The previous six-coil system provided only static field compensation with coil currents locked during an experimental session. As well as the large static components, significant temporal variation of the remnant magnetic field by as much as 10's of nT over the period of a single experiment (10-20 minutes) can occur depending on the environment. This drifting background field level has a detrimental effect on the open-loop OPM sensor gain. This gain varies by more than 10% over the full dynamic range of ±5 nT 13 . For accurate MEG measurements it is desirable that the sensor gain varies as little as possible, and as a result dynamic compensation of the background magnetic field is also required for sensor operation. Iivanainen and colleagues have described a field-nulling system which allows the operation of a QuSpin OPM array inside a poorly shielded environment by continually updating currents in a series of coils to account for dynamic variations in the remnant field 15 . In their work, dynamic nulling prevented significant sensor gain changes over the course of an experiment and reduced the peak-to-peak drift of the field from 1.3 to 0.4 nT 15 .
The three coils here were designed for use at a city-centre site located directly above a London Underground railway line. This environment produces unpredictable, temporally varying magnetic fields large enough to send OPMs outside of their operational range and frequent enough to make even a 10-minute scan unfeasible. To correct for this, a dynamic PID-based feedback controller was developed to continually update coil currents during a session in order to null reference measurements of B x , B y and B z . Empty room recordings with and without the dynamic nulling were compared whilst noting the number of trains which passed beneath the lab during these recordings.
Discrepancies between the predicted and measured coil performance may be caused by interactions between the coils and the high-permeability mu-metal walls of the MSR. This was not considered in the process of coil design here, though similar interactions have been incorporated into gradient coil design methods for MRI 16 . It has been shown that this interaction can be modelled by a method of images approach 17,18 . This was utilised to simulate the effects of interactions on the efficiency and homogeneity of the designed coils when operated inside a MSR. Results are compared to the measured performance and the impact this may have on future coil designs is discussed.

Coil Design Theory
We can characterise a current distribution J that is confined to an x-y plane using a two-dimensional stream function S(x, y), chosen such that ∇ × = z J S (since ∇.J = 0). The wire paths of a coil design that represents the current distribution correspond to the contours of the stream function. Planar current distributions which produce a uniform field along one of the Cartesian axes have shared symmetry with current distributions that generate specific field gradients, and this sharing of symmetry also applies to the stream functions. For example, in the case of a Helmholtz coil, two co-axial loops of radius a carrying the same current and placed a distance a apart in z will produce a uniform field B z along the z direction as shown in Fig. 2a. In this instance the stream function required to produce the field is symmetric in x and y, and also symmetric in z (as the loops carry the same current). If the current in one of the loops is reversed with respect to the other (so the stream function is now anti-symmetric in z, and symmetric in x and y) then the result is a gradient in B z along the z direction (i.e. a dB z / dz coil) as shown in Fig. 2b. If the two loops are driven independently the currents can be set to produce both a field, and a field gradient simultaneously as shown in Fig. 2c. This involves applying the sum of the currents needed to produce a field and a field gradient to one loop and the difference to the other loop.
Considering the symmetry of the coil stream functions required to generate each of the three uniform fields, the x-y symmetry is the same for B x (anti-symmetric in x, symmetric in y and anti-symmetric in z) and G x (anti-symmetric in x, symmetric in y and symmetric in z) coils, B y (symmetric in x, anti-symmetric in y and anti-symmetric in z) and G y (symmetric in x, anti-symmetric in y and symmetric in z) coils, and B z (symmetric in x, symmetric in y and symmetric in z) and G z (symmetric in x, symmetric in y and anti-symmetric in z) coils. Through a combined design process and careful control of the applied currents, three coil pairs can therefore be used simultaneously to generate the uniform fields and the associated gradients.
Here, we used a harmonic minimisation approach based on previously described expressions for planar magnetic field coils [19][20][21] to design the coil pairs. The two-dimensional Fourier transform of the magnetic field x y over the x-y plane at a position z in the region −a < z < a between the two planes can be written as: x y x y r r r x y k a 0 r where k r = (k x 2 + k y 2 ) 1/2 and  S k k ( , ) x y is the two-dimensional Fourier transform of the stream function. The upper/ lower terms define the z-symmetry where the stream function has the same (symmetric in z)/opposite (anti-symmetric in z) sign on each plane.
The stream function is then parameterised as a two-dimensional Fourier series confined to the region < x y L , ( = ± z a) on the two coil planes, so that 20 : where coefficients α n , β n , γ m and δ m are used to weight the different harmonics in the series. By applying symmetry conditions to the stream function, a sub-set of the harmonic combinations in Eq. (2) can be used to design each These stream function equations are written for ease of notation in the form Continuing with the case of the balanced B x /G x coil, the contribution to the field B x in the x-direction, b xj (r i ), from the j th component of the stream function can be expressed (defining the z-symmetry to be antisymmetric, i.e. we choose cosh in Eq. (1)) using Eq. (1) as where twoDFT denotes the two-dimensional Fourier transform. The contribution to the field from the gradient component G x in the x-direction, g xj (r i ), from the j th component of the stream function can also be expressed (defining the z-symmetry to be symmetric, i.e. we choose sinh in Eq. (1)) using Eq. (1) as i xj x j x y k a r i x y 0 , This can be substituted into Eqs (4) and (5) to find the field at position r i due to each component of the stream function. Similar calculations can be performed by imposing the symmetry conditions needed for the other coils. www.nature.com/scientificreports www.nature.com/scientificreports/ λ-coefficients are then chosen to minimise a functional containing three parts 20 : the spatial deviation of the field term (B x ) over an array of target points, the spatial deviation of the field gradient term (G x ) and the power dissipated in the coil. The functional is written as: Here, from the B x term B t (r i ) is the desired (or target) field at position r i and is the calculated field at r i . For the G x term g t z i is the desired (or target) field chosen to vary linearly with position z i and is the calculated field at r i . The set of position vectors r i = 1 to I define target points within the volume at which a homogeneous field or field gradient is required. The linear weighting α applied to the field gradient minimisation can be altered to balance the strength of the field and gradient terms. P is a tuneable power dissipation term which can be upweighted by increasing the weighting coefficient ω to reduce the complexity of the designed coils 20 .
Here the functional is minimised by choosing the weights λ j which satisfy dF/dλ j = 0. The set of derivatives can be cast as a set of linear simultaneous equations in matrix form, whose solution is found by identifying the pseudo-inverse matrix 7 . The wire paths of the coils are then extracted as contours of the optimised stream function.

Methods
Coil design and construction. MATLAB (The MathWorks Inc.) was used to produce designs for three coils allowing simultaneous compensation of B x and G x , B y and G y and B z and G z . The dimensions of the system were limited by the size and layout of the MSR which also contains a SQUID-based MEG system (CTF, Coquitlam, BC, Canada). Final bi-planar dimensions were chosen to be a = 0.7 m and L = 0.85 m as shown in Fig. 1. The number of harmonics, power term weighting and field/gradient term weightings were chosen manually to produce homogeneous fields and field gradients (to within ± 5%) over a central volume of 40 × 40 × 40 cm 3 . The magnetic field for each coil was evaluated over a 2.5-cm-resolution, 3D grid using the elemental Biot-Savart law. A numerical value of the field, or field gradient homogeneity, H was computed by assigning a value of 1 to every element in the 3D grid where the homogeneity condition was met, and a 0 where it was not. The homogeneity is then reported as a percentage of the number of 1's versus the total number of grid points. The balanced B x /G x and B y /G y , coils were designed using 36 harmonics with N = M = 6, and the field was evaluated over I = 320 target points. The balanced B z /G z coil was designed with 36 harmonics with N = M = 6 and the field was evaluated over I = 256 target points. Note that fewer target points are required to limit the field variation in the B z /G z coil due to the high degree of rotational symmetry of these coils about the z-axis.
Contours of the optimised stream function were connected to form continuous wire paths, as shown in Supplementary Information. The coils were mounted on two sheets of medium-density fibreboard (MDF) measuring 1.8 × 1.8 m 2 with legs which place the centre of the coils 1.1 m from floor level. Coil designs were printed and then hand-traced onto thin paper sheets, these sheets were attached to the board with wallpaper paste. Copper wire (enamelled, diameter of 0.56 mm) was laid on each path and fixed in place using masking tape. The total thickness of each coil layer is around 1 mm. The three coils were layered by repeating these steps.
Each planar coil unit was connected to a coil driver (4 V, low-noise, six drivers in total, one per coil unit), which was controlled using a 16-bit NI-9264 DAC voltage output module interfaced with LabVIEW (National Instruments (NI) Corporation, Austin, TX). Resistors were added in series to each coil unit so that fields of around 40 nT or field gradients of around 25 nT/m could be produced from a 4 V output.
Prior to their use in balanced mode, the coils were connected in standard, bi-planar arrangements (with equal or equal and opposite current in each plane 7 ) to test their performance when producing each of B x and G x , B y and G y and B z and G z . The coil efficiencies (field or field gradient per unit current) were measured by using an array of three Gen-1 QuSpin OPMs operating in a mode where the field is measured only along the sensor's radial axis. These OPMs were then oriented so as to measure the field component of interest.
A linearly increasing current was applied to each bi-planar coil and the change in field measured as a function of this applied current. For B x /G x and B y /G y the OPMs were positioned approximately at the centre of the bi-planar coils and spatially separated by 8 cm in the z direction (such that the field was measured at roughly z = −8 cm, 0 cm and +8 cm at x = 0 cm, y = 0 cm). For the B z /G z measurements, the long axes of the sensors were lined up in the measurement plane, which meant that the field had to be sampled at locations with a greater spatial separation in the z direction (field measured at z = −17 cm, z = 0 cm and z = + 17 cm at x = 0 cm, y = 0 cm). For comparison, the fields and gradients produced by the bi-planar coils were also calculated using the elemental Biot-Savart law. These measurements were carried out in the MSR at the University of Nottingham (Vacuumschmelze, Hanau, Germany) comprised of two layers of mu-metal and one layer of aluminium.

Comparing design methods.
To compare the performance of the balanced coils designed using the combined method to that of coils designed just to produce a uniform field or field gradient, three coils were investigated: a coil designed to produce B x ; a coil designed to produce G x ; a coil designed to produce B x and G x . For each case the number and location of target points and the weighting of the power dissipation term were kept constant. The dimensions of the coils were a = 0.7 m and L = 0.85 m as above. The field and gradient produced by these coils in a bi-planar arrangement was calculated using the elemental Biot-Savart law. Fields were then normalised to their values at the centre to allow evaluation of the deviation of the field and gradient over a 40 × 40 × 40 cm 3 volume (evaluated over a 2.5 cm resolution 3D grid). The homogeneity H is also reported for each case, using the measure described in section 3.1. As before, less than ±5% deviation over the full volume is desired. (2019) 9:14196 | https://doi.org/10.1038/s41598-019-50697-w www.nature.com/scientificreports www.nature.com/scientificreports/ Static field nulling. Although in their standard measurement mode QuSpin OPMs have an operational range of ± 5 nT they can also be operated in a "field-zeroing" mode 13, 14 . This mode produces the zero-field environment which is required for sensitive measurements based on changes in transparency of the vapour cell around a zero-field resonance 22 . Each OPM contains a set of three orthogonal "on-sensor" coils sited around the vapour cell. In field-zeroing mode currents are applied to the on-sensor coils under feedback control to maximise the optically-measured zero-field resonance. Once optimised currents have been found, the amplitude of the static magnetic fields in two directions perpendicular to the laser beam can be inferred from the applied current and the known field per unit current of each on-sensor coil. This field is reported in the range 1 to 50,000 pT.
Using the field-zeroing process four Gen-1 QuSpin OPMs were used to measure the three vector field components, at two positions spatially separated in the z-direction by ~40 cm, as shown in Fig. 3A. The measurements from the sensors are combined to form estimates of the magnitudes of B x , B y and B z and of G x , G y and G z (for example, G x is found through subtraction of two measurements of B x ).
To automatically null the field over the reference array, two proportional integral derivative (PID) controllers are driven in tandem and used to identify the optimal currents to pass through the coil windings in each plane. Taking the example of nulling the B x and G x fields , the reference array produces two measures of B x (b x1 and b x2 ) at locations which are spatially separated in z. The first controller varies a parameter ∆ in order to reduce b x1 − b x2 , which provides a measure of the field gradient G x . The second controller varies a parameter ∑ in order to reduce b x1 + b x2 , which provides a measure of B x . Currents I 1 and I 2 are then supplied to the two different coil planes, where I 1 = ∑ + ∆ and I 2 = ∑−∆. The balancing of the terms results in simultaneous reduction of the field and field gradient towards zero. A LabVIEW-based controller was developed to implement this algorithm by linking to the reference array and coil drivers. Once nulling is completed the PIDs are switched off, and stable currents are applied to the coils. This process is outlined in Fig. 3B. field mapping. To quantify the reduction in field and field gradient, field maps were taken over a single plane before and after the coils were used for field nulling. A 40 × 40 cm 2 region in the x-z plane lying between the reference sensors (which were separated by 40 cm in z) was mapped. By taking two measurements using the field-zeroing mode of a single Gen-1 QuSpin OPM at each position on a 10-cm grid (changing the orientation of the sensor to allow the field to be measured in three dimensions) over the plane, a representation of the static field can be obtained. The data were linearly interpolated (using the MATLAB interp2 function) onto a 1-cm resolution grid for visualisation, and field maps produced to show changes in magnitude and spatial variation.
The average field is calculated as the mean over each point in the interpolated grid with the associated standard deviation given as the error. A measure of the gradient is found from the difference in field (∆B) divided by the difference in z position (∆z). The difference of the field values measured at z = −20 cm and z = + 20 cm were taken so that ∆z = 40 cm. These measurements were carried out in the MSR at the University of Nottingham.
Dynamic field nulling. The balanced bi-planar coils described here were intended for use in an OPM-MEG system to be operated inside a MSR (Vacuumschmelze, Hanau, Germany) comprised of two layers of mu-metal and one layer of aluminium, which is sited in the basement of the Wellcome Centre for Human Neuroimaging at University College London (UCL). This city centre site is close to the Piccadilly Line of the London Underground as it runs between Russell Square and Holborn underground stations at a depth of around 33 m (maps generated from Transport for London data available here -https://www.dansilva.co.uk/down-underground). Given the MSR is in the basement of the Wellcome centre, and using Citymapper (www.citymapper.com) to estimate the distance to the line at ground level (around 85 m), we estimate the railway line to be ~90 m away from the MSR. As a consequence of proximity to the passing underground trains the background magnetic field in the MSR varies by several nT over a 10-minute period. This drift pushes the OPMs close to the edge of their dynamic range To drive a single coil two PID controllers are used. One uses Σ to drive the sum of the field at the two locations (driving for example the B x field down towards zero). The second PID uses ∆ to drive the difference of the field at the two locations (the gradient G x ) towards zero. Currents which are the sum and difference of these values are then applied to the planes. and causes gain changes. In order to address this issue, we developed a dynamic PID-based set-up that could be used to control the currents in the balanced, bi-planar coils to keep the ambient field within the sensors' operational range during recordings, minimising any impact on their performance.
An array of 15, Gen-2 QuSpin OPMs were placed in an empty 3D printed scanner-cast 11 . Three of these sensors, which were placed in slots oriented in the x, y and z directions so that they measured the field components B x , B y and B z along their radial axes, were selected as reference sensors. The sensors were operated in the single-axis measurement mode with the dynamic range set to ±5 nT and data were recorded using a National Instruments 16-bit NI-9205 ADC interfaced with a LabVIEW controller. The voltage outputs from LabVIEW for driving the coils were controlled using the NI-9264 DAC. The controller sampled data from all 15 sensors at a 6,000 Hz sampling rate and the data was saved in chunks of 50 samples. For each of these 50-sample chunks the data from the three reference sensors was selected and passed to a digital, first-order, low-pass Butterworth filter with a 1 Hz cut-off frequency (chosen such that the controller was addressing only the low frequency fluctuations in the remnant magnetic field). The average of these 50 samples was then taken.
The averaged values were then fed into three PID controllers (one for controlling each of the B x , B y and B z fields). These controllers were manually tuned by observing the system response to variation in the proportional and integral gains (derivative gains were set to zero) to provide stable nulling without interfering with each other, as cross-talk can be introduced between the coils and reference sensors. For example, the effect of poor coil or reference sensor positioning could lead to a measurable effect of the B x coil on the B y reference sensor, in turn leading to an over-correction on the B y coil which impacts on the B x reference sensor and so on. A well-tuned loop can readily compensate for these 'cross-talk' effects. The PID output voltages were applied to the coils on the two planes such that uniform field components were produced (i.e. voltages of the same sign were applied to the coils for generating B z , and of opposite sign for producing B x and B y ). Data from the remaining 12 sensors were also recorded. During each recording, data was taken from the Transport for London website for trains which ran on either of two routes: eastbound from Holborn to Russell Square and westbound from Russell Square to Holborn. Whilst this does not provide the exact time that trains came in closest proximity to the MSR, it does indicate the number of trains that crossed beneath the lab during the recordings.
A recording with the nulling off was taken between 17:22 and 17:32, while a recording with nulling-on was made between 17:35 and 17:45. These times fall during the 'rush hour' , during which there is heavy train activity (trains arrive roughly every 2 minutes on the Piccadilly line around this time), thus providing a challenging environment for the measurements.
Investigating interactions with the shielded room. The fields that are produced by the coils may interact with the walls of the MSR and result in deviations from the predicted performance. At a boundary between materials of different relative magnetic permeability μ r , the following conditions apply: to obey Ampere's law the tangential component of the magnetic field H must be continuous across the surface, and to obey Gauss's law the normal component of B must also be continuous 23 . Since B takes the form B = μ 0 μ r H inside the material and μ r can be >80,000 for mu-metal, while μ r ~ 1 in air, the field must abruptly change direction at an interface between air and mu-metal to satisfy the two boundary conditions 23 , and in particular on the air-side of the boundary the tangential components of the B and H fields must be approximately zero. Under these circumstances the field variation within the MSR can be evaluated by addition of a set of virtual mirror currents produced by simple reflection of current-carrying elements of the coils at each of wall of the room. For magnetostatics' problems, this "method of images" approach is not very sensitive to the exact value of μ r , just requiring that μ r » 1 since deviations from the field variation it predicts scale as (μ r ) −1 which is negligible for mu-metal 24,25 . To ensure that the boundary conditions are properly fulfilled, reflections have to be calculated recursively so that the mirrored currents are themselves mirrored. Here, the effect that the mu-metal walls have on the form of the field produced by the balanced coils operating inside the MSR has been modelled via the addition of a truncated set of mirror currents 17,18 .
Considering a single current element, reflections are made at each of the internal faces of the shielded room and the direction dl reversed according to = − .l l l n d d d 2( ) where n is the unit vector perpendicular to the plane in which the element was reflected. After consideration of each internal face of the MSR, this results in the production of 6 mirrored current elements. Each of these mirrored elements is then also mirrored by reflection at the five walls of the MSR that were not involved in the original reflection. This process should be repeated ad infinitum fully to satisfy the boundary conditions on the magnetic fields at the MSR's internal walls. Accounting for redundancies there are 6 first order mirror elements, 18 second order mirrors, 38 third order mirrors and so on. The additional field produced by each of the mirrored elements is calculated using the elemental Biot-Savart law. Here we only considered reflections up to 3 rd order. This order was chosen as the change in coil efficiency (at the centre of the planes) with reflections up to second order was found to be around 30%, while the additional contribution from third order mirrors only modified the efficiency by a further ~5%.
We conducted a simple experiment to test the mirroring method in the UCL MSR, which has internal dimensions of 3.85 × 2.4 × 3 m 3 . A single coil plane was placed parallel, and as close as possible, to one of the 3.85 × 2.4 m 2 walls of the MSR. The coil plane was offset from the centre of the MSR in by 0.6 m and −0.2 m in x and y and due to the legs of the wooden support structure the centre of the coil was displaced by 25 cm in the z-direction from the proximal wall of the MSR. A 0.1 Hz sinusoidally oscillating voltage with amplitude of 1 V (equivalent to 8.1 mA current) was applied to the B z /G z coil on this plane and the z-component of the magnetic field was measured at a position 5 cm displaced from the centre of the coil in the z-direction away from the proximal wall of the MSR (but centred in x and y) using a single-axis fluxgate magnetometer (Fluxmaster, Stefan Mayer Instruments, Dinslaken, Germany). The coil plane and fluxgate magnetometer were then further displaced from the wall in steps of 20 cm in the z-direction until the setup could not be moved any further due to the presence of a cryogenic MEG system (CTF, Coquitlam, BC, Canada). A final displacement of 85 cm from the wall was achieved. At each position the fluxgate magnetometer was similarly positioned relative to the coil. For comparison the Scientific RepoRtS | (2019) 9:14196 | https://doi.org/10.1038/s41598-019-50697-w www.nature.com/scientificreports www.nature.com/scientificreports/ z-component of the field was calculated, with, and without the effect of mirroring by applying the elemental Biot-Savart law to the B z /G z coil and its reflections.
The effects of the interactions on the bi-planar coils was then investigated by simulation. The field was calculated over a 40 × 40 × 40 cm 3 , 2.5 cm resolution 3D grid at the centre of the two coil planes and the results compared to measurements of the coil field or field gradient per unit current. Three cases were studied: first, when no mirrors were considered, second when the coils were positioned centrally in the MSR at the University of Nottingham and mirrors up to third order were included, and finally when the coils were positioned as they are during an experimental session and mirrors up to third order were considered. For each case a box plot of the full range of fields calculated at all grid positions was produced to show the effects of the mirror currents on both the homogeneity and the mean level of the produced field.
The effects were investigated for the three coil pairs when producing B x , B y and B z . The MSR considered here has internal dimensions of 3.85 × 2.4 × 3 m 3 and during an experiment the centre of the bi-planar coil set is displaced from the centre of the room by + 0.2 m in y and −0.25 m in x (with zero offset in z).

Results
Coil designs. Figure 4 shows wire paths (i) for the (a) B x and G x , (b) B y and G y and (c) B z and G z coils respectively. The elemental Biot-Savart law was applied to the wire-paths to evaluate the spatial variation of the field (ii) or field gradient (iii) produced, when the coils are wound in a bi-planar arrangement. The field was calculated over a 3D array of points covering a region of 40 × 40 × 40 cm 3 . Contours are shown over a plane at z = 0 m; |x|, |y| < 0.2 m for the case of B x , B y , B z , G x , and G y and over the plane at z = 0.025 m; |x|, |y| < 0.2 m for the G z case. The fields were normalised to the value at the centre of the two planes (or a value positioned slightly off-centre when the field is zero at the centre). The homogeneity of all the produced fields was H = 100%.  Table 1 shows calculated and measured values for the bi-planar coil efficiencies and single plane inductances. The errors quoted for the field generating coils are the standard deviation from the mean measured across three OPMs spatially separated in z, while the errors in the gradient generating coils are the differences in the two gradient values calculated from the three OPM measurements. The inductances were calculated from the stream functions using previously published expressions 26 . Comparing design methods. Figure 5 shows for coils optimised to produce (a) B x , (b) G x and (c) and both B x and G x (i) designs for coils, (ii) contours of B x and (iii) contours of G x . The field and gradient variation is  www.nature.com/scientificreports www.nature.com/scientificreports/ calculated using the elemental Biot-Savart law and is plotted over a plane at z = 0.20 m; |x|, |y| < 0.20 m. This plane is deliberately chosen to be off-centre to highlight the decrease in homogeneity with distance from the centre of the region between the two coil planes.
Comparing the homogeneity of the non-optimised field (e.g. the G x variation produced by the B x coil) shows the combined method improves homogeneity, but with an increase in design complexity. The B x coil had a G x homogeneity of H = 76% and the G x coil had B x homogeneity of H = 78% compared to H = 100% for the optimised field components. The largest deviations of the field and gradient at z = 0.20 m in Fig. 5a(ii-iii) are +2% and −13% respectively. The field and gradient in Fig. 5b(ii-iii) deviate by −9% and +2% respectively. The field and gradient in Fig. 5c(ii-iii) deviate by +3% and −4% respectively. Only the combined design satisfies the condition that both the field and field gradient should be homogenous to within ±5% over a central 40 × 40 × 40 cm 3 cubic region.
field mapping. Figure 6 shows maps of the spatial field variation in (a) B x , (b) B y and (c) B z inside the MSR before (i) and after (ii) the automated nulling process was applied, along with a bar chart (iii) showing the reduction in average field strength. The mean magnitude of the field vector was reduced by a factor of 26 from 23 nT to 0.9 nT. The reduction in the measured individual field components with associated standard deviation from the mean is shown in Table 2. The maps of the field variation are not as smooth as are expected to be found in a source-free region, as a result of inaccuracies in manually setting the sensor position and orientation, and the limited accuracy of measurement that can be achieved when operating the QuSpin sensors in field-zeroing mode.
Dynamic field nulling. Figure 7a,b show the time-series data from the 15 sensors placed in an empty scanner-cast inside the empty MSR. Figure 7a shows data recorded over a 10-minute period with the dynamic nulling off, while Fig. 7b shows data recorded with dynamic nulling on. Vertical lines show the approximate times that a train arrived at either Russell Square or Holborn stations during the two different recording periods, thick lines show two trains arriving at the same time. During the nulling-off recording, 10 trains were reported in the 10-minute period, while for the nulling-on recording 7 trains were reported.  Figure 6. Quantifying the coil performance by measuring the field in a single plane. (a-c) show for B x , B y and B z respectively: (i) a map of the field before nulling is applied; (ii) a map of the field after nulling is applied; iii) the average field magnitude with error-bars showing the spatial standard deviation of measurement before (red) and after (blue) nulling is applied.
The maximum peak-to-peak variation of the time-series data across the 15 sensors was reduced from 2.37 nT to 0.13 nT by the dynamic nulling. The maximum standard deviation from the mean of the sensor outputs across sensors was reduced from 0.501 nT to 0.012 nT. Figure 7c shows the median power spectral density estimate across all 15 sensors. To compute this estimate, the data were split into three non-overlapping sections of 200-s duration. For each section, flattop windowing (using the MATLAB function flattopwin) was applied prior to calculation of the power spectral density via fast Fourier transform. The results were then averaged over the three sections for each sensor and the median taken across all 15 sensors. No additional processing (e.g. forming synthetic gradiometers or linear regression) was applied to the raw magnetometer traces to remove other sources of interference. As a result, the noise spectra are mostly dominated by interference sources, which exceed the noise level specifications of the OPM sensors (https://quspin.com/products-qzfm/).  Table 2. Comparing the measured magnitudes of each of the field components before and after field nulling. Errors are standard deviation from the mean.  Figure 7d shows the shielding factor of the system computed as the ratio of the PSDs shown in Fig. 7c for the nulling-off and nulling-on cases. The shielding factor is reported in decibels (dB). Figure 8 shows the measured B z -field amplitude as a function of the z-displacement of the coil plane from the wall of the MSR, errorbars are calculated as the standard deviation of the amplitude over ten cycles. The expected enhancement of the measured field amplitude with proximity to the wall of the MSR is seen. Good agreement with the calculated data is also evident, indicating that the mirroring method provides a reasonable model of the interaction between the coils and the walls of the MSR. Figure 9 shows box plots representing the mean and spatial variation over a central 40 × 40 × 40 cm 3 volume of the calculated (A) B x , (B) B y and (C) B z fields per unit current (i.e. coil efficiencies) for three cases, where: (i) no mirror currents were considered; (ii) the bi-planar coils were assumed to be centrally positioned within the MSR and mirror currents up to third order were considered; (iii) the bi-planar coils were positioned off-centre, as they would be during an experiment, and mirror currents up to third order were considered. The measured experimental value of the efficiency and its associated standard deviation is shown as a red star with error bar for each case. As indicated previously, the experimental values of B x B y and B z were averaged over three OPMs spatially separated in z. If the experimental value lies within the range of the box plots this is a good indication of the model being broadly accurate. High spatial homogeneity is desirable for producing a controlled field over the subject.

Investigating interactions with the shielded room.
For the case of the B x coil with no mirroring of currents, the median calculated field over the grid was 0.72 nT/mA with a range of 0.05 nT/mA (maximum value 0.7 nT/mA and minimum value 0.75 nT/mA). Positioning  www.nature.com/scientificreports www.nature.com/scientificreports/ the coils at the centre of the MSR and adding mirror currents reduced the median calculated field to 0.64 nT/ mA with a range of 0.06 nT/mA (max/min 0.61/0.67 nT/mA). Realistic positioning of the coils within the MSR and accounting for current mirroring further reduced the median efficiency to 0.62 nT/mA with a range of 0.21 nT/mA (max/min 0.67/0.46 nT/mA). The measured value of the efficiency of 0.57 ± 0.01 nT/mA lies within the increased range of field values. The simulations show that the field homogeneity decreases when the effect of the MSR is included. The maximum deviation of the field over the 40 × 40 × 40 cm 3 volume that was considered increases from 5% with no consideration of mirror currents to 26% when the effects of the mirror currents are included for the realistically positioned coils.
Similar results were obtained for the B y coil with the median simulated field decreasing from 0.66 nT/mA with a range of 0.05 nT/mA (max/min 0.64/0.69 nT/mA) when mirror currents were not considered to 0.54 nT/ mA with a range of 0.17 nT/mA (max/min 0.6/0.43 nT/mA) when the coils were realistically positioned and the effects of the mirror currents were included. This range includes the measured value of 0.47 ± 0.02 nT/mA. The field inhomogeneity is increased from 5% to 22% for this coil when the mirror currents are included In the case of the B z coil, the median field increased from 3.90 nT/mA with a range of 0.11 nT/mA (max/ min 3.97/3.86 nT/mA) to 4.27 nT/mA with a range of 0.61 nT/mA (max/min 4.45/3.84 nT/mA) on inclusion of mirror currents with realistic coil positioning. The upper bound of the error on the experimentally measured value of 3.60 ± 0.30 nT/mA lies within this range, but the general trend of an increased efficiency is not observed. The inhomogeneity is increased from 5% to 10% for this coil when the effect of the mirror currents is taken into account.

Discussion
Coil design and construction. By introducing a design and control method to produce coils which provide field nulling of six field terms using only three coil designs, construction of a field nulling system for use in wearable MEG is considerably simplified. The performance of the new field-nulling system is comparable to that obtained previously using a six-coil system 7 .
Although simpler designs such as the 'Helmholtz-cage' system described by Iivanainen and colleagues do not require complex wirepaths or use of large biplanar coils the subject has to remain fixed within a cubic enclosure 15 . A compromise has to be made upon system design between the freedom of motion and the quality of the MEG data which is obtained. Shielded rooms with larger remnant fields and gradients and larger temporally varying fields will pose a greater challenge for field nulling and as a result a fixed system may be simpler to implement. Additional compensation of higher order terms as well as improved dynamic nulling could allow for subject movement in poorer shielded environments.
To fully compensate for the uniform field and field gradients 8 terms (3 uniform fields and 5 gradients) are required, as the divergence and curl of the field are both zero in the current-free region between the planes. The system described here only produces 6 field terms. In this bi-planar arrangement the additional dB x /dy and dB x / dx terms have no shared symmetry with the uniform coils (or with one another) and would have to be individually constructed. Designs for these coils have been detailed previously, but they were not constructed 7 .
Although the number of coils that have to be wound is halved with this approach, the process of manual construction is still a cumbersome task and will likely result in errors such as inaccurate wire placement, misalignment from the wallpapering of the coils and warping of the wood. Although the PID approach allows for some of these errors to be taken into account, and the field nulling performs well, 3D printing, printed circuit board etching or water cutting techniques could potentially be used to increase the accuracy of coil winding and decrease labour requirements. With these techniques, the complexity of wire paths that could be produced would also increase.
Static and dynamic field nulling. Though the system described here provides nulling which allows for movements by a seated subject, additional coils could also be constructed to provide better nulling and potentially allow a greater range of movements. These could include coils to generate the remaining two field gradients (dB x /dx and dB x /dy), individual, higher order spherical harmonics or more complex local patterns of field variation. We showed that dynamic nulling of uniform field components could be implemented using the balanced bi-planar coils, and this allowed the OPMs to be operated continually in a magnetically noisy environment, meaning that experiments could be completed without saturation of sensor outputs due to interference from a nearby underground railway line.
The relatively low mutual inductance (M) and self-inductance (L) of the coils (all <400 μH) and the low frequency of operation (<120 Hz) meant that the effects of inductance were not important when driving the coils since MdI/dt and LdI/dt were both always much smaller than IR (where I is the current and R is the coil resistance). The slight increase in noise level which is introduced in the nulling-on case (see Fig. 7c,d) produces a negative shielding factor at some frequencies above 1 Hz. This noise results from the limited dynamic range of the field nulling electronics used here. To compensate for the remnant static field components the field nulling system has to be able to produce tens of nT; the NI-9264 DAC module has a ± 10 V dynamic range and a 16-bit resolution, the least significant bit is therefore equivalent to 20/2 16 ~ 0.3 mV. In an example system, where 40 nT is required from 10 V this is equivalent to a field of roughly 1 pT. Rounding errors in setting the voltage in dynamic nulling consequently lead to additional field variation at the pT level at the sensor. In the low frequency range (<1 Hz) the environmental field variation is much larger than 1 pT and the dynamic nulling performs well, producing a significant shielding factor. However, at higher frequencies the environmental field noise is lower in magnitude and the additional field noise from the limited dynamic range of the nulling system becomes the dominating factor. This additional noise is apparent in the presented data, but as the effects are systematic and picked up by all sensors, synthetic gradiometry could be used to attenuate them. www.nature.com/scientificreports www.nature.com/scientificreports/ To attenuate these effects in future field nulling systems a 24-bit DAC could be used to increase the dynamic range. The development of a two-circuit driver, one optimised for static nulling (producing e.g. 40 nT from 10 V) and one for dynamic nulling (producing e.g. 4 nT from 10 V) with the sum of the currents applied to the coils could be developed.
The speed of the LabVIEW based controller is fundamentally limited by the loop rate of the software. More advanced controllers, such as field programmable gate array microcontrollers or analogue solutions will be required for increased nulling performance at higher frequencies.
Although here the reference sensors were mounted in a scanner-cast during the MEG recording they would have to be placed behind the subject as shown in Fig. 3 in a real experiment, since with head-mounted reference sensors subject movements would reorient the reference sensors away from the x, y and z directions. This setup would allow the dynamic nulling to be readily expanded to include compensation of the field gradients G x , G y and G z through measurement of the full vector field on both sides of the head. This has not yet been implemented due to the limited number of sensors available for use in the reference array. A more complex arrangement of the reference array and feedback mechanism will most likely be required to provide an optimal dynamic null over an entire head array of OPMs. Flexible arrays, body mounted reference sensors and prior knowledge of the static magnetic field at different points combined with optical tracking could be used to provide a stable null as the subject moves within the MSR. Such improvements in the performance of coil technology for field compensation will be needed to facilitate new avenues of OPM-MEG research, such as experiments involving spatial navigation and direct social interaction between individuals.
Interactions with the shielded room. Our brief investigation of the use of mirror currents to model interactions with the MSR revealed a change in median efficiency and reduction in homogeneity of the produced field which broadly agreed with the experimental measurements (Figs 8 and 9). The results imply the need to incorporate the effects of the high-permeability mu-metal into future coil designs, but it should be noted that the PID approach to both static and dynamic field nulling will seek to find the best solution without relying on a priori information.
Comparison of the boxplots in Fig. 9 shows that, for all three coils, the effect of the MSR on the field homogeneity is much less when the coils are positioned at the centre of the room. The off-centre positioning of the coil results in an asymmetric distribution of the additional field components from the mirror currents which make the field inhomogeneity larger than the ±5% condition targeted at the design stage. This asymmetry will also affect the efficiency and homogeneity of the gradient components, distorting the spatial distribution of the produced field. Unfortunately, the coils cannot be positioned centrally in the MSRs at UCL and Nottingham, as a ~500 kg cryogenic MEG system is sited at the centre of each room. Increased displacement of the coils from the centre of the MSR in the x and y directions would further effect the spatial distribution of the mirrored current elements; the closer a coil plane is to the wall, the greater the asymmetry of the resulting field distribution.
The experimental values were taken over single lines with a 16 cm spread for the B x and B y cases and a 34 cm spread for the B z case instead of being measured over the entire 40 × 40 × 40 cm 3 region. This may account for some of the discrepancy between experimental and theoretical measurements. In general, the interaction between the coils and the mu-metal walls should be more thoroughly investigated in order to incorporate the effect of the MSR when designing future coils. A more thorough experimental test of the range of validity of the simple mirror approach for static and dynamically-varying fields would be worthwhile, complemented by electromagnetic calculations using finite element analysis software.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on request.