On-scalp magnetocorticography with optically pumped magnetometers: Simulated performance in resolving simultaneous sources

Currently, the gold standard for high-resolution mapping of cortical electrophysiological activity is invasive electrocorticography (ECoG), a procedure that carries with it the risk of serious morbidity and mortality. Due to these risks, the use of ECoG is largely limited to pre-surgical mapping in intractable epilepsy. Nevertheless, many seminal studies in neuroscience have utilized ECoG to explore domains such as visual perception, attention, auditory processing, and sensorimotor behavior. Studies such as these, occurring in patients with epilepsy rather than healthy controls, may lack generalizability, and are limited by the placement of the electrode arrays over the presumed seizure focus. This manuscript explores the use of optically pumped magnetometers (OPMs) to create a non-invasive alternative to ECoG, which we refer to as magnetocorticography. Because prior ECoG studies reveal that most cognitive processes are driven by multiple, simultaneous independent neuronal assemblies, we characterize the ability of a theoretical 56-channel dense OPM array to resolve simultaneous independent sources, and compare it to currently available SQUID devices, as well as OPM arrays with inter-sensor spacings more typical of other systems in development. Our evaluation of this theoretical system assesses many potential sources of error, including errors of sensor calibration and position. In addition, we investigate the influence of geometrical and anatomical factors on array performance. Our simulations reveal the potential of high-density, on-scalp OPM MEG devices to localize electrophysiological brain responses at unprecedented resolution for a non-invasive device.

Introduction
Although other invasive electrophysiological approaches such as microarrays (Truccolo et al., 2011) and stereotactic EEG exist (Nagahama et al., 2018), the gold standard for high-resolution mapping of the cortex is invasive electrocorticography (ECoG). Because of the high risks of morbidity and mortality (Rolston et al., 2015), ECoG is generally only used in situations such as pre-surgical planning for epilepsy. Despite this, many seminal studies have utilized ECoG in patients with epilepsy to also investigate basic neurocognitive processes. Due to the prevalence of temporal lobe epilepsy, the auditory and language systems have been studied extensively using ECoG. A study of language processing revealed that even a relatively simple word recognition task activated seven areas encompassed by standard subdural grids, with up to four areas simultaneously active (Zheng et al., 2021). While one ECoG study revealed distinct neuronal populations performing high frequency amplitude and low frequency phase tracking of auditory speech (Zion Golumbic et al., 2013), another identified distinct areas for expressive and receptive speech (Towle et al., 2008). ECoG has also been used to investigate other domains. A study of visual processing revealed a complex interplay between phase and frequency across cortical sites (Miller et al., 2010). During imagined movement, ECoG revealed distinct neuronal populations generating alpha and beta rhythms in sensorimotor cortex (Stolk et al., 2019). Without ECoG, research has shown that even the simplest appearing processes in the cortex may consist of many simultaneous subprocesses; for example, the N100 auditory evoked potential consists of at least six underlying simultaneous components (Naatanen and Picton, 1987). Studies investigating correlation structure of traditional and microarray ECoG have found that independent signals can be detected at sub-millimeter resolution (Maier et al., 2010;Bullock and McClune, 1989), and that correlations between electrodes are largely due to the spatial extent of local field potentials (Rogers et al., 2019). Taken together, these studies suggest that ECoG can resolve multiple, closely spaced, simultaneous independent sources.
While MEG is a unique modality for measuring electrophysiological activity of the brain with both high temporal resolution (greater than 100 Hz), and moderate (5-10 mm) spatial resolution, whole-head MEG systems are not capable of distinguishing among closely spaced independent sources of similar amplitude. To understand why that is the case, factors affecting the spatial resolution of MEG systems need to be considered. Traditional SQUID based MEG systems require a fixed helmet with substantial thermal insulation, which inevitably creates a gap between the scalp and the sensor and imposes a fixed inter-sensor distance. Besides sensor array geometry, other factors affecting the spatial resolution are the inherent sensor noise, the accuracy of the forward model, the signal-to-noise ratio (SNR) of the measurements, the presence of interfering sources, and the cortical geometry (including cortical surface orientation) apart from any co-registration issues.
Technological progress in recent years have led to the development of magnetic field sensors called optically pumped magnetometers (OPMs), which have the potential to overcome the limitations of the current SQUID based MEG devices. OPMs measure magnetic field changes by measuring the transmission of laser light through a gas of alkali metal atoms (such as Rubidium, 87 Rb). Fluctuations in the magnetic field affect the transparency of the vapor, which results in changes in the amount of light at the photodetector site. OPMs can be placed directly on the scalp, thereby decreasing the distance between the sensors and the cortical sources. This arrangement allows sampling of higher spatial frequencies of cortical activity, dramatically improving SNR, and thus improving spatial resolution (Boto et al., 2016;Borna et al., 2020). In addition, decreasing the distance between adjacent sensors can also improve the ability of MEG to resolve and accurately discriminate among multiple independent sources. OPM systems have additional advantages including the lack of need for cryogens and the ability to construct wearable systems that can move with the participant (Boto et al., 2018). Most work developing OPM MEG arrays has focused on whole-head systems to replace traditional SQUID MEG. In contrast, a system designed to accurately discriminate among multiple, closely spaced simultaneous sources would need to evaluate activity for a restricted region of the cortex, but at greater spatial resolution (analogous to ECoG grids). While a whole head array at this density may not be currently feasible, OPM arrays sampling a portion of the cortex at high density are presently achievable and are thus the focus of this investigation.
Prior authors have simulated various aspects of OPM array performance and their comparison to commercial SQUID MEG arrays. The simplest simulation involves measuring signal to noise; simulating a single dipole predicts a ten-fold increase in SNR for an on-scalp whole head system compared to a traditional axial gradiometer SQUID MEG system (Iivanainen et al., 2017). Correlation of the simulated time course for a source with its reconstructed time course is another convenient measure of source reconstruction accuracy. For cortical sources, correlations between simulated and reconstructed sources were greater than 0.9 everywhere on the cortical surface for an OPM array, with lower correlations observed for an axial gradiometer SQUID system (Boto et al., 2016). Additionally, the simulated OPM system exhibited less signal leakage into neighboring dipole locations (Boto et al., 2016). Another metric involves correlating source topographies to determine localization accuracy. Using this metric, average peak position errors were similar between radial OPM magnetometers (5.1 mm), SQUID magnetometers (7.9 mm) and SQUID planar gradiometers (5.2 mm) (Iivanainen et al., 2017). Another OPM simulation study showed that errors in the position of sensors (i.e., 2-6 mm) result in large (7.5-20%) errors in the reconstructed topography of the lead fields, i.e., fields generated by each source at the location of each sensor. Likewise, errors in the orientation of the sensor sensitive axis (5-15°) also result in large (4.8-14.7%) errors in the reconstructed topography (Zetter et al., 2018). Another straightforward metric is to measure the distance between the location of a simulated dipole source and the peak of the linearly constrained minimum variance (LCMV) beamformer source reconstruction. One study calculated this distance to be 7.9 mm for an array of OPM magnetometers, while mean errors of 10.7 mm and 11.4 mm were calculated for SQUID magnetometers or a combined magnetometer/ planar gradiometer SQUID system, respectively (Zetter et al., 2018). Fewer simulations have investigated the performance of systems using multiple simultaneous sources. In the presence of five interference sources, the correlation between reconstructed and simulated time courses for a source of interest was found to be 0.7 for an OPM system and around 0.5 for an axial gradiometer SQUID system (Boto et al., 2016). These same simulations showed that the improvement in correlation relies strongly on accurate lead field estimation (Boto et al., 2016). Another simulation in the same report showed that an OPM system could differentiate two dipoles with orthogonal time series that were 2.5 mm apart, while a SQUID system could only resolve sources 6.5 mm apart (Boto et al., 2016). Importantly, while this report did utilize multiple simultaneous sources, the authors did not investigation how the discrimination of sources was affected by sensor array characteristics or calibration errors. In interpreting these studies, it should be noted that metrics related to source localization accuracy of beamformers are dependent on signal to noise and may vary between recovered sources in an image.
The purpose of this study is the evaluation of a proposed OPM sensor array designed specifically for magnetocorticography (MCoG), and its comparison with existing magnetoencephalography (MEG) instruments, with a focus on the discrimination of simultaneous sources. This simulation study differs from and extends previous work in several ways. First, our simulations are based on a spatially dense OPM sensor array that, like ECoG, is optimized for resolving activity within a small patch of cortex, while most previous work has focused on whole head OPM arrays with more widely spaced sensors. One study has simulated high density whole head OPM arrays but does not comprehensively address errors inherent in OPM systems (Boto et al., 2016). Furthermore, whole brain arrays with high resolution are not currently feasible; one study determined that 500 sensors would be required to spatially discriminate sources at 6 mm (Tierney et al., 2020). Nevertheless, given that we focus this investigation on the cortical area directly under the array, a dense whole-head array would be expected to have similar performance for cortical sources. The second novel contribution of this study is that we focus on differentiating simultaneous independent sources instead of focusing on single sources, as has been the focus of most prior studies. We will quantify performance in part by defining an independent source resolution metric (ISR) as a measure of the ability of a sensor array to discriminate source activity at biologically plausible intervals. Third, unlike the work of previous investigators, we recognize that a crucial factor limiting the resolution of independent sources is the presence of interfering source activity. While prior investigators have considered single or few (Boto et al., 2016;Brookes et al., 2021) interfering sources, we have developed a brain noise model to which sensor noise and test sources can be superposed. This model is validated by comparison with task-free brain activity measured with existing MEG sensor arrays. Lastly, we assess the effect of factors involved in the design of a sensor array, including the gap between the array and the scalp, and errors in the known positions and orientations of the sensors, as well as errors attributable to phenomena such as crosstalk and perturbation of the true sensor axis.

Methods
The purpose of the simulations is to characterize the performance of a densely packed array of OPMs that can be developed using existing technology. Because commercially available MEG systems currently in use utilize traditional low-temperature SQUID systems, we include the simulated performance of a SQUID MEG system for comparison purposes. Thus, simulations are based on SQUID and OPM sensors currently used in the National Institutes of Mental Health MEG Core Facility: a CTF 275 whole head MEG system (CTF MEG, Coquitlam, B.C., Canada) and FieldLine V2 OPM sensors (FieldLine, Inc., Boulder, CO). While other commercially available sensors exist, potentially with the advantage of multiple axes of sensitivity, we focus on a sensor offering closed loop mode on a single axis to avoid nonlinearity and reductions in sensitivity inherent in a multi-axis sensor. Because we are not attempting to fully characterize the SQUID system, we evaluate it "as built" rather than assessing performance in the presence of errors.

Noise model
The ability of a magnetometer sensor array to resolve multiple independent sources is limited by the presence of interfering neural sources (the brain noise), sensor noise, sensor geometry, and cortical geometry. In this section we describe the methods used for noise simulation. Sensor noise for both SQUID and OPM sensors was simulated using a Gaussian random sample generator with zero mean and standard deviation equal to the specified root mean square (RMS) noise. We assumed a noise density of 5 fT/√Hz for the SQUID sensors and 10 fT/√Hz for the OPM sensors. Brain noise was generated by an empirical model. We present our final model and brief description of our validation here and present the full verification of the model in the Supplemental Materials. All simulations were performed using in-house code developed in C (distributed as part of the SAM suite https://github.com/ nih-megcore/SAMsrcV5), with python used for visualization. Realistic brain noise was simulated by 2000 dipoles randomly located within a hemispherical shell ranging from 60 to 80 mm in radius. Because most MEG signals are generated by cortical pyramidal cells, we constrained sources to a spherical shell rather than the full volume. The 20 mm shell thickness was chosen to capture the depth below the skull of most sources, given gyrification of the cortex. The amplitude of noise sources ranged from Q min = 1.0 pA-m/√Hz to Q max = 1.0 nA-m/√Hz, according to the following power law, where p = 6. Q i = (i/N) p Q max − Q min + Q min All brain noise dipole sources had a random orientation in the plane tangential to the surface of the sphere. In addition, we assume coupling among brain noise sources. This is modeled by applying a modulus parameter of M = 10 such that every group of 10 consecutive dipoles had the same time series (scaled by the previously defined random noise amplitude). In order to verify that our model approximated ongoing brain activity, we placed a single simulated test source in a source localized resting state dataset, as well as in the noise shell of our simulated geometry. We compared noise density, width of the beamformer reconstruction of the test source, and the eigenvalue spectra from real and simulated data.

Geometry and simulation characteristics
The simulation geometry is shown in Fig. 1. The origin was established as the center of the hemispherical shell containing the simulated brain noise sources. For the OPM system ( Fig. 1A), the OPM sensors are arranged in a rectangular array 7 sensors wide and 8 sensors long to achieve a roughly square array, given the 13 × 15 mm footprint of the commercial OPM sensors. As a baseline, all simulations use a 13 × 15mm spacing of the sensors, although as detailed later, we explore the effects of varying this density. Published values for scalp to cortex distance are approximately 15 mm for the motor cortex (Lu et al., 2019). Adding extra distance for a cap or other device holding the array, as well as to compensate for the standoff of the 87 Rb cell within the fixture containing them (on the order of 4-5 mm), we set the OPM array at a radius of 102 mm (22 mm above the brain noise shell). For simulations using the CTF-275 MEG system (Fig. 1B), we used sensor positions derived from our system at NIH, which has 272 functional sensors. The maximum z-coordinate for the CTF array was approximately 110 mm. This would place the apex of the CTF array approximately 30 mm above the noise shell. Note that this would represent a best-case-scenario where the head was as high as possible into the MEG helmet; the majority of the sensors would be substantially further away from the sources.
The array of test source dipoles were randomly positioned within a 50 × 50 mm square, projected onto a spherical surface within or on the outer surface of the noise shell. As introduced in the last section, brain noise sources occupied a hemispherical shell with radius ranging between 60 and 80 mm. The simulated test dipoles were assigned random orientations in a plane parallel to a plane tangential to the sphere.
The duration of all simulations was 300 s. The sample rate for the OPM sensor simulations was 175 Hz, twice the bandwidth limit of the sensors (approximately 350 Hz), while a typical sampling frequency of 600 Hz was used for the CTF simulations. Each test dipole was assigned a unique random time series, scaled by a source magnitude of 25 nA-m RMS (similar to the 22 nA-m sources utilized in (Zetter et al., 2018)), unless otherwise stated. The magnetic field observed at each sensor was simulated by summation of brain noise, sensor noise, and the test dipole sources. The field generated by dipole sources was calculated assuming a spherical homogeneously conducting space using the Sarvas equation (Sarvas, 1987). Note that because we are assuming a sphere of activity, only tangentially directed source currents contribute to the magnetic field external to the spherical conductor.

Array calibration errors
In order to understand the performance of both idealized OPM arrays and realistic OPM arrays, we simulated various potential errors that may occur at the sensor level. These include crosstalk, angular errors to the sensitive axis, and gain error, each of which we will explain in turn. OPM sensors make use of three orthogonal onboard field nulling coil pairs to establish zero field in the 87 Rb cell, to keep the sensor in a domain of relative linearity. When operated in closed loop, the field produced by these nulling coils is dynamically altered to maintain zero field throughout a recording. While closed-loop operation will ensure a linear response, the field produced by a nulling coil on one sensor may also influence the field at neighboring sensors. This phenomenon is known as crosstalk, and it is important to understand its effects on data quality. While crosstalk cannot be eliminated, it can be mitigated in part using self-shielding nulling coils (Nardelli et al., 2019) and calibrating measured data using the inverse of the measured crosstalk matrix. We examined the effects of crosstalk ranging from 1 to 10% on the measured signal. In this simulation, sensors at the minimum distance on the array experience the full percentage effect of crosstalk, while the effect is reduced according to the inverse of the distance cubed for sensors at further distances.
In addition to crosstalk arising from closed-loop operation, there is also modulation crosstalk. The sensing axis of each OPM is established by applying high-frequency modulation to the z-axis Helmholtz coils in the sensor and demodulating that signal to measure the opacity of the 87 Rb vapor to the probe laser. Crosstalk between the z-axis coils of adjacent OPM sensors due to this modulation can perturb the sensing axis, contributing to inaccuracies in the measured signal. We explored the effects of orientation error, which we also refer to as angular error, from 1 to 4°, by giving each sensor a fixed angular offset from the ideal axis of the sensor, at a random angle in the tangential plane. While this is also a phenomenon due to crosstalk, it is static throughout the measurement provided that the sensors do not move relative to one another Sensors may also differ in their precise gain values. We modeled errors in sensor gain by scaling the observed signal by the sensor by a random percentage up to the maximum gain error specified. We explored the effects of gain errors up to 10%.

Perturbations to array and simulation geometry
The second group of potential issues affecting reconstructed signals involve the geometry of the array. While ideally, the OPM sensors would rest on the surface of the scalp, this may be impractical due to concerns with sensor temperature (the vapor cell is maintained at an internal temperature of approximately 150°C; see (Shah et al., 2007)). While custom 3D printed head casts would conform to a participant's head, a fixed helmet which can be used on many different participants may result in gaps between the scalp and the face of the sensor. We simulated the effect of scalp-sensor gap up to 12 mm to account for this scenario. Flexible caps conforming to the participants head can be used as an alternative to head casts and helmets; however, this may result in inaccuracies in the precise localization of sensors (Hill et al., 2020). Finally, we examined the spacing of sensors. For our default sensor spacing, we assumed sensor housings of 15 × 13mm (same as FieldLine OPM sensors), and that sensors were placed immediately adjacent to one another. To investigate the effect of using densely packed arrays vs. arrays resembling whole head coverage, we simulated array spacings ranging from 10 × 8mm to 20 × 18mm (including 13×11, 14×12, 15×13, 16×14, and 17 × 15mm). Although the current footprint of the FieldLine v2 sensors is 15 × 13mm, we extrapolated to smaller spacings to more comprehensively show the dependence of SNR and ISR on sensor spacing.
In addition to these intrinsic OPM characteristics, we also examined other factors that affect the recorded signal, namely the number and magnitude of test sources and the cortical depth at which they are located. To explore the relative performance of OPM compared to SQUID MEG systems in resolving independent sources, we located the test dipoles in a 50 mm square area and varied their number from 2 to 128, and their magnitude from 5 nA-m to 30 nA-m. Given that our noise model utilized a 20 mm thick shell of noise sources, we can model the effects of cortical depth by varying the location of the test sources from the outer radius to the inner radius. While the thickness of the grey matter ribbon is considerably less than 20 mm, gyral folding should result in deeper sources of interest. The response of the array to sources at varying cortical depth will be in part due to the increased distance between sources and sensors, decreasing the spatial frequencies captured. In addition, the response will be influenced by the noise sources between the test dipoles and the sensors.

Source reconstruction
For each iteration of array and dipole arrangement, the reconstructed time series for each test source was obtained using a linearly constrained minimum variance (LCMV) scalar beamformer. The beamformer requires the calculation of the covariance matrix of the (simulated) measurements, C = m(t)m(t) T , and the lead-field matrix. The lead field matrix L r is the field generated by a unit dipole at position r in the two orthogonal directions tangential to the sphere. The elements b r are the forward solution vector components for each dipole location. (Note that for a realistic head model, the matrix L would have three columns corresponding to the cardinal directions x, y and z).
The scalar beamformer proceeds in two steps. First, the dipole source orientation is estimated by solving the generalized eigensystem of the ratio of projected source power to projected noise power.
L T C −2 Le k = λ k L T C −1 Le k k = 1, 2 In the preceding equation, λ k and e k are the eigenvalues and eigenvectors, respectively. Next, the forward solution for the source is computed for the projection of the eigenvector corresponding to the maximum SNR (largest eigenvalue) with the lead field: b r = L r e max Lastly, the scalar beamformer coefficients giving maximum SNR are: Thus, the estimated source time series and projected noise become: S r (t) = w r T m(t) and N r = w r T σ 2 W r respectively, where σ 2 is the sensor noise variance. In addition to beamformer estimates of the source time series at the known location of each test source, we computed images from the signal-to-noise ratio (source power divided by the projected noise), at each 1 mm voxel in a 60 × 60 mm grid according to the following equation: SNR r = w r T Cw r w r T σ 2 W r All software was written in-house as part of a larger suite of imaging analysis tools.

Test metrics
Several metrics were used to evaluate the simulated arrays. In order to provide a measure of how well different sensor arrays can discriminate between simultaneous sources, we calculated a metric we refer to as Independent Source Resolution (ISR). For each test source, we computed the correlation between the simulated and reconstructed time courses, which is a measure of the fidelity of reconstruction of the sources. In addition, to measure the effects of signal leakage, we calculated the correlation between the simulated time course of one test source with the reconstructed time courses of all other test sources. Thus, a correlation matrix can be formed, with each element being the correlation between the simulated time course (along one axis of the matrix) and the reconstructed time course (along the other axis of the matrix). The ISR is then defined as the mean on-diagonal correlation minus the mean off-diagonal correlation. The use of correlations between reconstructed and true test sources has been used previously to investigate the reconstruction of a source given brain inference (Boto et al., 2016). Our particular implementation is novel, as we utilized a much larger number of simultaneous sources, allowing richer characterization of array performance.
The ISR calculation, as well as the correlations used to derive the ISR are presented in several ways. First, for each test array we calculated the mean ISR and SNR across all test sources for 25 different source configurations (by incrementing the initialization seed of the random number generator). Next, for a single configuration of test sources we plotted the full correlation matrix, where, in an ideal scenario, each of the simulated dipole sources would be reproduced with high fidelity with no leakage of signal (and thus the correlation matrix would resemble an identity matrix). In addition, the correlation values between test dipole pairs were plotted against the physical distance between them. Correlation will be reduced as test sources are further apart, and the degree of decay of the correlation with distance can be used to infer how well individual sources can be resolved.
While these metrics are continuous, it may be useful to define an acceptable threshold for source discriminability. Boto et al. (Boto et al., 2016) suggest that any two sources displaying shared variance under 50%, corresponding to a correlation coefficient of 1/√2 = 0.7, are discriminable. Using the concept of shared variance, we would define a more stringent threshold, and propose that any two sources sharing less than 10% variance can be differentiated, corresponding to a correlation coefficient of 0.32 (i.e indicating acceptable leakage). In addition, we would suggest that an accurate reconstruction of a source should show at least 80% shared variance with the true source, corresponding to a correlation coefficient of at least .89 (indicating acceptable reconstruction accuracy). Thus, a desirable ISR would be greater than 0.57, although it remains useful to consider reconstruction accuracy and leakage separately.
Because our time courses for each test dipole are reconstructed at the a priori known location of the test source, a source whose reconstruction is shifted in space will exhibit potentially artifactually lower SNR. To assess how different array characteristics affect the accuracy of source localization, we also calculated and plotted two-dimensional images on the spherical surface containing the test sources (again for a single configuration of test sources). As with the reconstructed time courses, these images were constructed using an LCMV beamformer, and we plot projected source power divided by projected sensor noise power (SNR).
Finally, it should be noted that in a real system, multiple inaccuracies will be present simultaneously, further degrading performance, although we would expect these errors to add linearly rather than to interact. In addition, as the SNR of signals increases, the beamformer will more closely approximate a point source, leading to reduced output SNR if the measurement vector and forward model do not agree (Vrba, 2002). We would expect this phenomenon to also be independent of other influences on ISR and SNR.

Baseline performance of the CTF and OPM array
In Fig. 2, we show simulation results for the CTF MEG array and the theoretical OPM array under ideal conditions -no crosstalk or other sensor calibration errors, superficial sources, and minimum scalp to cortex distance. The OPM array clearly outperforms the CTF array in the resolution of 16 independent 25 nA-m (ISR OPM = 0.91, ISR CTF = 0.63). Correlation matrices clearly show few off-diagonal correlations, indicating less signal leakage for the OPM array. The enhanced resolution of sources is also evident in the beamformer images, as is the increased signal to noise, with a maximal SNR of 65.81 in the OPM image compared to a maximum SNR of 7.39 in the CTF image.

Errors of sensor calibration: crosstalk, sensitive axis error, and gain error
The simulated effects of crosstalk are shown in Fig. 3. Given that crosstalk is one of the most widely discussed issues with OPM sensors, we expected it to have a significant effect on all parameters we measured. In contrast, we found that crosstalk related to closed-loop operation only had a relatively minimal impact on ISR and the ability to differentiate simultaneous independent sources (Fig. 3 B and C). As is evident in Fig. 3 A, ISR for the OPM array remained well above that found for the CTF array, although SNR did decrease markedly as a function of crosstalk magnitude. In the reconstructed source images calculated using the simulated data ( Fig. 3 D), all sources were visible and properly localized at crosstalk magnitudes of up to ten percent.
In contrast to crosstalk induced errors resulting from closed-loop operation, angular error in the sensitive axis of the sensor, introduced by modulation crosstalk, had a profound effect for the proposed OPM array on all examined metrics, as shown in Fig. 4. Mean ISR was reduced to unacceptable levels between 2 and 3 degrees of error (ISR OPM for angular error of 3° = 0.51), while SNR levels of the OPM array were reduced below the SNR of the CTF array at around 3 degrees of error (See Fig. 4 A). The primary driver of the reduction in ISR was due to a sharp decrease in the correlation between each source's measured and simulated time courses (Fig. 4 B and C). These reductions were due in part to a mis-localization of the sources due to the axis error. The ISR is calculated from the recovered signal at the location of each test dipole; if the test dipole position is shifted due to errors in the geometry of the sensors, the ISR would appear more reduced than if the recovered source time course from the shifted location was utilized. As is evident in Fig. 4 D, even a one degree perturbation in the sensitive axis results in shifts of the peak locations of the recovered sources.
Similarly, errors in gain (shown in Fig. 5) resulted in dramatic reductions in both ISR and SNR. ISR for the OPM array was reduced to levels simulated for the CTF array between 4 and 5% gain error, while SNR was reduced dramatically to below that simulated for the CTF array between 3 and 4% gain error (See Fig. 5 A). As with axis errors, rather than inducing signal leakage between independent sources, errors in gain reduced the fidelity of the reconstructed time series (see Fig. 5 B and C). Gain errors of 5% resulted in images without distinct sources visible (Fig. 5 D).

Effects of sensor and array geometry: scalp-to-array gap and sensor position errors
Even for gaps of up to 12 mm, scalp-to-array gap had relatively low impact on the ability to accurately reconstruct source time series (See Fig. 6 A-C; for a 12 mm gap ISR = 0.71).
Although ISR remained relatively high, the array with a large gap showed greater signal leakage, leading to a reduction in spatial resolution. Examining the reconstructed source maps (Fig. 6 D), it was readily apparent that sources appeared much less compact, and more diffuse, like the spread in activation observed for the CTF system. Nevertheless, both ISR and SNR remained above that observed for the CTF system. It should be noted that this is in part a demonstration of the on-scalp OPM array to capture the high spatial frequencies present in our test dipole array.
Sensor position errors additionally had a profound effect on the ability of the theoretical array to discriminate independent sources (Fig. 7). Similar to errors in the sensitive axis of the sensors, this is partially due to mis-localization of sources. Even deviations of 0.5 mm of the location of the sensors was enough to severely degrade the reconstructed signals and render images of the test dipoles nearly indistinguishable from noise sources (Fig. 7 D).
In order to fully demonstrate the effects of sensor spacing, we evaluated the results using both 16 and 32 simulated sources. As the spacing between sensors increased, the profile of the spatial reconstruction of the sensors broadened. When 32 test sources were evaluated, mean SNR increased as sensors were spread further apart due to the overlap of source reconstructions (Fig. 8). However, when 16 test sources were evaluated, mean SNR at source locations decreased due to mislocalization of source peaks. ISR remained relatively high for all configurations of the OPM array, dropping to that of the CTF system only for 32 sources and the maximum sensor spacing (20 × 18mm).

Performance of OPM and SQUID systems: number, magnitude, and depth of test dipoles
The final set of parameters we assessed were aspects of the test dipoles themselves: their number, magnitude, and depth. A comparison of the CTF array and the theoretical OPM array at distinguishing multiple sources is shown in Fig. 9. The theoretical OPM array was capable of accurately discriminating among 32 dipolar sources (ISR OPM = 0.75). In addition, although the ISR was below our threshold for 64 sources (ISR OPM = 0.53), examining the plot of correlation vs. distance between test dipoles shows that sources at separations of at least 10 mm could be reasonably distinguished. In contrast, the CTF system struggled to resolve more than 16 sources (For N = 32, ISR CTF = 0.43).
When test dipole magnitude was reduced, the OPM system was able to distinguish 16 sources at the lowest magnitude tested (5 n-Am; ISR OPM = 0.59, see Fig. 10). In contrast, the CTF system had difficulty distinguishing sources and accurately localizing them at source strengths below 15 n-Am (10 n-Am; ISR CTF = 0.55). When the test dipole strength was increased, the ISR observed for the theoretical OPM system approached ISR = 1 for 16 independent test dipole sources, with SNR increasing as well. In contrast, for the CTF system, the maximal ISR of approximately 0.62 was achieved for a 20 nA-m source, and further increases in source strength did not increase ISR. Likewise, SNR peaked close to 5 for a source strength of 15 nA-m, then slowly decreased as source strength increased. Further examination showed that reconstructed source power increased with dipole magnitude, but the projected noise power in the area of each source also increased and became more focal. It is unclear if this is a true property of the system or a product of our simulation parameters and geometry.
Finally, we examined the effect of cortical depth of sources in the noise shell (Fig. 11). While cortical depth had a profound effect on both the ISR and SNR of the theoretical OPM array, the ISR value was greater for the OPM array than the CTF array for all cortical depths (although within the margin of error for the deepest sources, Fig. 11 A). SNR remained higher for the OPM array at all depths, although deeper sources were poorly localized ( Fig.  11 A and D). For the OPM array, cortical depth decreased the correlation between each test source and its reconstructed time series and increased the correlation with the remaining test dipoles. In contrast, for the CTF array, while the within-dipole correlation also decreased, the correlation between sources was relatively unaffected (Fig. 11 B and C).

Discussion
In this manuscript, we evaluated an OPM array designed specifically for magnetocorticography and compared its performance in discriminating multiple independent sources with an existing SQUID MEG instrument. We assessed the effects sensor calibration issues, including crosstalk, angular errors in the sensitive axis of the sensors, and errors in sensor gain. We also evaluated errors in the presumed location of the sensors, the effect of the gap between the sensors and the scalp, and the effect of sensor spacing. Finally, we assessed the performance of our array in discriminating sources when we altered the number, strength, and depth of test dipoles.
Overall, we found that crosstalk due to closed-loop operation affected source reconstruction less than we expected. Even with 10% crosstalk for nearest neighbor sensors, the OPM array had superior ISR and SNR compared to the CTF system. This is an important finding, given that crosstalk will be exacerbated as sensors are placed closer to one another, which is required for high-resolution applications. In Boto et al. (Boto et al., 2018), crosstalk was measured in an array of sensors 3-5 cm apart and was on the order of 3%. The sensors developed by FieldLine incorporate self-shielded zeroing coils, which reduces crosstalk from 8% to 0.5% for sensors that are immediately adjacent to one another (Nardelli et al., 2019).
In contrast, the OPM system was highly sensitive to angular errors in the sensing axis of the sensor, as well as errors in gain introduced in part due to modulation crosstalk. Angular errors of 3% or gain errors of 5% resulted in unacceptable reductions in ISR. Similarly, we showed that the positions of sensors must be known to a high degree of accuracy, with errors no larger than 0.25 mm, consistent with prior published results (Zetter et al., 2018;Iivanainen et al., 2019). One caveat to this is that it may be difficult to calibrate and localize sensors with this degree of accuracy, given inherent noise in the system and factors such as crosstalk and gain errors (Duque-Munoz et al., 2019), as well as the influence of ambient transverse fields (Iivanainen et al., 2019). While authors have experimented with data-driven methods to determine orientation offsets, only angular errors above 4° can be rejected in the best case (Duque-Munoz et al., 2019). These results indicate that flexible sensor array designs (similar to EEG caps) where sensors may move relative to each other during a recording, are likely unsuitable as modulation crosstalk would be expected to change.
Provided that the position of sensors is known, however, there is some flexibility in the sensor-to-scalp distance and the spacing of sensors. While minimum gaps between the sensors and the scalp, and between adjacent sensors results in the highest separability of simultaneous sources, gaps of 15-18 mm between the 87 Rb cells and sensor-to-scalp gaps of up to about 8 mm offer acceptable resolution, provided SNR is sufficient given the location and strength of sources. This is advantageous given that space for heat mitigation strategies such as airflow around the sensors may be required to for long duration measurements.
Given its higher signal to noise ratio, the OPM system outperforms the CTF SQUID system for test dipoles of every evaluated magnitude, although the difference in performance is most evident for the stronger dipoles. This is partially a consequence of the reduced sensitivity of OPM sensors (10 vs. 5 fT/√Hz). Similarly, examining depth of test sources within the noise shell, the OPM system does not substantially outperform the CTF system for the deepest sources in terms of ISR, although SNR of the OPM system remains greater than that observed for the CTF system. This is consistent with prior simulation work, showing that the advantage of OPM systems compared to SQUID systems becomes smaller as the depth in the brain increases (Boto et al., 2019). This is in part a consequence of the greater head coverage of the CTF array. We also examined the number of simultaneous sources. Conceptualizing brain activity during a task as a finite number of discrete sources is an oversimplification, but is a useful model, particularly for sensory tasks. As the number of sources increases, the reconstruction problem becomes underdetermined; our 56 channel array cannot mathematically provide a unique solution in the presence of more than 56 sources. The OPM array nevertheless performs quite well in comparison to the CTF array, which cannot differentiate 32 sources (ISR OPM = 0.75, ISR CTF = 0.43) The broader ramifications of the enhanced ability of OPM systems to discriminate simultaneous sources is readily apparent. Functional magnetic resonance imaging (fMRI), with the advent of ultra-high field scanners, has enabled the mapping of retinotopic (Olman et al., 2010), somatotopic (Martuzzi et al., 2014), and tonotopic (Norman-Haignere et al., 2013) representations. While somatotopic finger mapping has been performed using MEG, alternating stimulation between digits is required for precise localization, so that only one source is presumed active at a time (Jamali and Ross, 2013). Likewise, retinotopic maps, for instance of stimulus orientation, can also be reconstructed with conventional MEG, although with greater difficulty. One study, for instance, mapped three ranges of stimulus eccentricities, although only in V1, and mapping of the upper visual field and sulcal sources was unsuccessful (Nasiotis et al., 2017). Hierarchical visual processing in the visual system necessarily results in multiple simultaneous sources that may be possible to visualize using a dense on-scalp OPM array. Tonotopic maps have also been investigated using MEG (Falet et al., 2021), although subjects were highly heterogeneous and not all participants exhibited a gradient in receptive field from low to high frequency tones. In addition to sensory mapping, ultra-high field MRI has enabled laminar-specific imaging (Huber et al., 2018). There is evidence that the frequency of cortical oscillations may be layer specific (Sun and Dan, 2009), suggesting that the laminar specificity of MEG signals could be inferred indirectly. While Bayesian modeling has been used to differentiate MEG signals arising from deep or superficial cortical layers (Bonaiuto et al., 2018;Pinotsis and Miller, 2020), true layerspecific MEG imaging remains unattainable. Work using dense OPM arrays, however, could potentially enable a new frontier in the study of cortical lamina. In addition to the ECoG studies discussed in the introduction, other studies using EcoG revealed distinctive activation patterns for motor movements (Aoki et al., 1999), and successfully decoded hand movement and force (Flint et al., 2020), suggesting a clear application for a dense OPM array in brain computer interface. Likewise, EcoG recordings have been used to decode human speech, which may suggest another application for OPM arrays (Pasley et al., 2012;Kellis et al., 2010). Of course, a magnetocorticography array with similar resolution to EcoG has an obvious application to characterizing the seizure onset zone in epilepsy. Small studies have performed assessment of interictal activity using OPM arrays (Mellor et al., 2021;Vivekananda et al., 2020). Given the demonstration of these comparatively sparse arrays to localize spike activity, our dense array would be expected to perform quite well, at least for superficial cortical spikes. Commonly used ECoG arrays have an electrode spacing of approximately 10 mm, with spatial resolution probably in the millimeter range (Wilson et al., 2006). Depth electrodes, frequently utilized in presurgical planning, are thought to measure local field potentials within a sphere of radius 0.5-3 mm surrounding the electrode tip (Logothetis, 2003). Although spatial resolution is not straightforward to measure, given that it depends on both source number (Fig. 9) and source strength (Fig. 10), it is likely that our system would achieve millimeter resolution. Nevertheless, we do not explicitly model a comparison between ECoG and our OPM array, so definitive conclusions cannot be drawn.
There are several important limitations to this work. Regarding the sensor array, the relatively small size of the array restricts the field of view, and sources outside that field of view will add noise to the measurements. Although the most obvious solution would be to sample the entire brain with a high density array, this is not currently feasible in practice, although this option has been addressed in simulation studies (Boto et al., 2016). In addition to sources outside the field of view there is also aliasing from higher spatial frequency information, although our dense array would suffer less from this problem than arrays with larger sensor spacings (Tierney et al., 2020). As for simulation methods, we used a spherical simulation rather than a realistic head model. While this is an approximation, it simplified our geometry and ensured that our results were not confounded by individual cortical folding patterns. Our use of a finite number of equal magnitude sources is also an approximation of actual brain behavior but may nonetheless be relevant especially in the case of epileptic spike generation as well as simple sensory tasks. Our empirically derived model of brain noise may not accurately reflect underlying generators of ongoing brain activity; nevertheless, our inclusion of brain noise represents a strength of this analysis. We should also note that we specifically simulate the behavior of FieldLine sensors, which have a single axis of sensitivity and closed-loop operation for improved linearity. There are currently commercially available biaxial and triaxial OPM sensors, but these do not operate in closed-loop mode, making nonlinear response a potential limitation. In addition, increasing to two or three sensitive axes results in a reduction of the overall sensitivity of the sensor. Thus, we confined our investigation to single axis sensors, although it may be of interest to repeat these simulations for multi-axis sensors. Recent studies showed that the inclusion of triaxial sensors in an array enhanced uniformity of coverage (Boto et al., 2022) and improved the rejection of external interference and movement artifacts, although increasing the density of single axis OPMs produced better rejection of internal sources (Brookes et al., 2021). Triaxial sensors operating in closed loop with high sensitivity may represent a significant improvement over the currently simulated array, once such sensors are commercially available. Likewise, hybrid arrays of OPMs and EEG electrodes have also been proposed in the literature (Beltrachini et al., 2021), although we do not explore such configurations here.
In this manuscript, we demonstrated the ability of a theoretical OPM magnetocorticography array to resolve multiple independent sources to a degree currently achievable only using invasive ECoG. We simulated most effects which could be expected to degrade performance of an OPM array and demonstrated that OPM arrays remain superior to conventional SQUID MEG for a range of potential calibration errors or other experimental issues. These results should serve to convince the reader that OPM sensors permit the design of arrays that can achieve unprecedented resolution in a non-invasively, potentially transforming the field of human neuroscience.

Citation diversity statement
As proposed in (Zurn et al., 2020), we include a statement of gender diversity in our reference list. The gender balance of papers cited within this work was quantified using a manual gender determination from authors' publicly available pronouns. From the 42 references included in this work, the distribution of genders for the first and last authors were as follows: 71.4-76.2% were MM (N = 30-32), 7.1% were MW (N = 3), 14.3-16.7% were WM (N = 6-8) and 2.4% were WW (N = 3). The gender of the first authors of 2 Illustration of the simulation geometry for the A) magnetocorticography array simulation and the B) CTF-275 simulation. Left panes illustrate the geometry in 3D while the right panes illustrate a projection onto the y-z plane. Scattered green dots represent the noise sources, while test sources and sensors are depicted in red and blue, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2.
Comparison of the CTF and idealized OPM system with no calibration errors. Metrics shown are A) correlation between source time series plotted as a function of the distance between sources, B) correlation between source time series shown in matrix form, and C) LCMV beamformer image of the test source area. True test source locations are marked in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3.
Effects of crosstalk on evaluation metrics. A) mean ISR and SNR for varying degrees of crosstalk across 25 different random sets of 16 test dipoles for crosstalk values between 0 and 10%. B) Correlation between source time series plotted as a function of the distance between sources, C) correlation between source time series shown in matrix form, and D) LCMV beamformer images of the test source area for 1% and 10% crosstalk. Effects of angular error in the OPM sensitive access on evaluation metrics. Figure legend for individual panels as in Fig. 3; mean ISR and SNR is shown for angular error between 0 and 5°, other evaluation metrics (correlation [C] and LCMV beamformer images [D]) are shown for angular error of 1, 2, and 3°. Effects of gain error on evaluation metrics. Figure legend for individual panels as in Fig. 3; mean ISR and SNR is shown for gin errors between 0 and 10%; other metrics (correlation [C] and LCMV beamformer images [D]) are shown for 1, 2, and 5% gain error. Effects of an additional gap between the OPM sensor and the test sources. Figure legend for individual panels as in Fig. 3; mean ISR and SNR is shown for additional gaps between 0 and 12 mm; other metrics (correlation [C] and LCMV beamformer images [D]) are shown for gaps of 2 mm and 10 mm.  Effects of spacing between individual OPM sensors. A) mean ISR and SNR for different sensor spacings across 25 different random sets of both 16 (blue data points) and 32 test dipoles (green data points). Because commercially available sensors have rectangular rather than square faces, spacings are not symmetric, and range from 10 mm × 8 mm (not currently possible with any commercially available sensor) to 20 mm × 18 mm. B) Correlation between source time series plotted as a function of the distance between sources, C) correlation between source time series shown in matrix form, and D) LCMV beamformer images of the test source area are shown for spacings 10 mm × 8 mm, 15 mm × 13 mm, and 20 mm × 18 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 9.
Effects of the number of tests sources. A) mean ISR and SNR across 25 different random sets of 4-128 test dipoles, shown for both CTF (orange markers) and the theoretical OPM system (blue markers). B) Correlation between source time series plotted as a function of the distance between sources, C) correlation between source time series shown in matrix form, and D) LCMV beamformer images of the test source area are shown for both the CTF and theoretical OPM systems, for 8, 16, 32, and 64 test dipoles. Effects of the strength of the test dipoles on evaluation metrics. A) mean ISR and SNR across 25 different random sets of test dipoles ranging in magnitude from 5 nA-m to 30 nAm, shown for both CTF (orange markers) and the theoretical OPM system (blue markers). B) Correlation between source time series plotted as a function of the distance between sources, C) correlation between source time series shown in matrix form, and D) LCMV beamformer images of the test source area are shown for both the CTF and theoretical OPM systems, for 5 and 30 nA-m dipole strengths. Effects of the depth of the test dipoles inside the hemispherical noise shell, representing a simulation of cortical depth. A) mean ISR and SNR across 25 different random sets of test dipoles ranging in depth from 0 mm to 20 mm, shown for both CTF (orange markers) and the theoretical OPM system (blue markers). B) Correlation between source time series plotted as a function of the distance between sources, C) correlation between source time series shown in matrix form, and D) LCMV beamformer images of the test source area are shown for both the CTF and theoretical OPM systems, for cortical depths of 4 mm and 16 mm.


s2fieldsofstudy: Biology Medicine

text:
PET/CT Radiomic Features: A Potential Biomarker for EGFR Mutation Status and Survival Outcome Prediction in NSCLC Patients Treated With TKIs

Backgrounds Epidermal growth factor receptor (EGFR) mutation profiles play a vital role in treatment strategy decisions for non–small cell lung cancer (NSCLC). The purpose of this study was to evaluate the predictive efficacy of baseline 18F-FDG PET/CT-based radiomics analysis for EGFR mutation status, mutation site, and the survival benefit of targeted therapy. Methods A sum of 313 NSCLC patients with pre-treatment 18F-FDG PET/CT scans and genetic mutations detection were retrospectively studied. Clinical and PET metabolic parameters were incorporated into independent predictors of determining mutation status and mutation site. The dataset was randomly allocated into the training and the validation sets in a 7:3 ratio. Three-dimensional (3D) radiomics features were extracted from each PET- and CT-volume of interests (VOI) singularly, and then a radiomics signature (RS) associated with EGFR mutation profiles is built by feature selection. Three different prediction models based on support vector machine (SVM), decision tree (DT), and random forest (RF) classifiers were established. Furthermore, nomograms for estimation of overall survival (OS) and progression-free survival (PFS) were established by integrating PET/CT radiomics score (Rad-score), metabolic parameters, and clinical factors. Predictive performance was assessed by the receiver operating characteristic (ROC) analysis and the calibration curve analysis. The decision curve analysis (DCA) was applied to estimate and compare the clinical usefulness of nomograms. Results Three hundred thirteen NSCLC patients were classified into a training set (n=218) and a validation set (n=95). Multivariate analysis demonstrated that SUVmax and sex were independent indicators of EGFR mutation status and mutation site. Eight CT-derived RS, six PET-derived RS, and two clinical factors were retained to develop integrated models, which exhibited excellent ability to distinguish between EGFR wild type (EGFR-WT), EGFR 19 mutation type (EGFR-19-MT), and EGFR 21 mutation type (EGFR-21-MT). The SVM model outperformed the RF model and the DT model, yielding training area under the curves (AUC) of EGFR-WT, EGFR-19-WT, and EGFR-21-WT, with 0.881, 0.851, and 0.849, respectively, and validation AUCs of 0.926, 0.805 and 0.859, respectively. For prediction of OS, the integrated nomogram is superior to the clinical nomogram and the radiomics nomogram, with C-indexes of 0.80 in the training set and 0.83 in the validation set, respectively. Conclusions The PET/CT-based radiomics analysis might provide a novel approach to predict EGFR mutation status and mutation site in NSCLC patients and could serve as useful predictors for the patients’ survival outcome of targeted therapy in clinical practice.

INTRODUCTION
Lung cancer is the leading cause of cancer-related death worldwide. Each year, approximately 1.6 million people die of lung cancer, and its five-year survival rate ranges from 4% to 17% (1). Histologically, non-small cell lung cancer (NSCLC) is the most frequent pathological subtype, which accounts for about 85% of the cases. Although early-stage lung cancer patients have a higher postoperative survival rate, treatments of advanced NSCLC show a relatively low response rate and significant toxicity (2). With the advance of precision medicine and personalized treatments, targeted therapy of NSCLC plays an increasingly important role as a rising star and was demonstrated to effectively improve the survival prognosis of lung adenocarcinoma patients with EGFR gene mutations (3). A series of previous studies (4,5) have shown that patients with EGFR mutations exhibited longer overall survival (OS) and progression-free survival (PFS) than those with EGFR-WT when receiving tyrosine kinase inhibitors (TKIs) therapies. Additionally, regarding the most common sensitive mutations include exon 19 deletion (19DEL) and exon 21L858R, previous studies have demonstrated that patients with 19DEL mutations may have a greater survival benefit after TKIs treatment than those with 21L858R missing mutations (6,7). Therefore, NSCLC therapies underwent an innovative transformation when it was realized that the mutant status of epidermal growth factor receptor (EGFR) directly affected the effectiveness of EGFR TKIs. It is critical to identify the molecular profiling of EGFR status in advanced NSCLC prior to individualized targeted therapy.
At present, clinical gene mutation detection usually uses tissue or cytological specimens, which has some disadvantages, such as trauma, difficulty in sampling, high cost, and unavoidable temporal and spatial heterogeneity of tumors (8). Analysis of circulating cell-free tumor DNA (ctDNA) is considered to be another emerging method for assessing EGFR mutation status (9). However, studies have shown that the ctDNA test has a relatively high false negative rate in clinical application, and the price is relatively high (10,11). Therefore, there is an urgent need to develop noninvasive, simple, rapid, and reliable methods for gene mutation detection.
Radiomics is an emerging field in which a large number of quantitative imaging features are extracted from medical images to identify those most closely related to clinicopathologic, molecular, and genetic characteristics with the purpose of improving the diagnostic and prognostic accuracy (12). Although a series of works (13)(14)(15) have been reported to explore the potential relation between EGFR mutation status and radiomic features derived from CT images, only a few studies using PET/CT have been reported in this field. In the molecular imaging, it is often based on visual analysis or conventional parameters, maximum standardized uptake value (SUVmax), e.g., resulted in unideal predictive performance. Nevertheless, there is a lack of related researches integrating radiomics features with conventional semantic features. Moreover, previous studies mainly focused on the differentiation between EGFR-WT and EGFR-MT without involving the identification of specific mutation sites (EGFR-19-MT or EGFR-21-MT).
Therefore, the purpose of this study was to investigate whether radiomics features extracted from the same volume of interest (VOI) of PET and CT images combined with metabolic indexes and clinicopathological parameters could be used to predict EGFR mutation profiles and mutation site based on a triclassification method. Furthermore, we intended to predict survival benefits of NSCLC patients treated with TKIs.

Patient Selection
This study was approved by the institutional review committee of Harbin Medical University Cancer Hospital. Given the retrospective nature of the study design and the anonymity of patient information, the informed consent requirement was waived. A total of 313 histologically proven NSCLC patients were retrospectively enrolled who underwent pretreatment 18 F-FDG PET-CT scans in our hospital between January 2013 and June 2018. Inclusion criteria were as follows (1): pathologically confirmed NSCLC (2); PET-CT scans performed within one month prior to surgery or biopsy (3); no history of any antitumor therapy before scanning (4); no history of other malignancies (5); a single lesion with a maximum diameter ≥ 1 cm. Exclusion criteria were as follows (1): no genetic test for EGFR or unavailability of genetic test results (2) none or low FDG metabolism of pure ground-glass nodules (3) incomplete clinical data (4) difficulty in tumor margin delineation. Clinicpathological information was obtained through clinical medical record retrieval, including age, gender, pathological stage, location, adenocarcinoma predominant subtype, carcinoembryonic antigen (CEA), smoking history and tumor size. Metabolic data including SUVmax, mean standardized uptake value (SUVmean) and total lesion glycolysis (TLG) were also recorded. The dataset was randomly assigned in a 7:3 ratio to the training cohort and validation cohort. Study design and patient allocation are shown in Figure 1. All cases in the training cohort were used to train the classification model, while cases in the validation cohorts were used to independently evaluate the model's performance.

EGFR Mutation Detection
Specific gene mutation information is confirmed by performing genetic testing on tumor tissue samples obtained by surgical resection or biopsy by an experienced physician. The mutation sites of four exons (exon [18][19][20][21] in the coding region of the EGFR gene were detected by real-time PCR. If any exon mutation was identified, the tumor was classified as EGFR-MT, otherwise considered as EGFR-WT.

Image Acquisition
All patients fasted for more than 6 hours before scanning, and were tested blood glucose levels, which were kept below 11.0 mmol/L. The image acquisition was performed using the discovery VCT 64 PET/CT system (GE Healthcare, Milwaukee, USA). A 3.78 MBq/kg dose of FDG was administered intravenously. Approximate one hour later, whole-body CT scanning was performed with a standardized protocol consisting of 120 kV, 140 mA, and 3.75 mm slice thickness. Then, for PET, the images acquisition time was 2 minutes per bed position. Image reconstructions were performed based on the 3D ordered subset expectation-maximization algorithm (2 iterations and 17 subsets).

Image Analysis, Tumor Segmentation and Radiomics Feature Extraction
The PET/CT images were analyzed by two radiologists blinded to the clinical and pathological results, (Reader 1, M.W and Reader 2, M.P with 15-and 20-years' experience in the interpretation of PET/CT images, respectively). The metabolic parameters were measured by drawing a region-of-interest (ROI) on the axial PET image based on a threshold of 40% of SUVmax using commercial  software (PET VCAR; GE Healthcare, USA). Any disagreement was resolved by consensus. SUVmax was defined at the highest value on one pixel with the highest counts within the ROI (16). The overview of radiomics workflow is displayed in Figure 2. Axial PET and CT digital imaging and communications in medicine images obtained from the picture archiving and communication system were applied for tumor segmentation. The tumor lesion was delineated separately on axial PET and CT images using LIFEx software (open-source software; www. lifexsoft.org/index.php). All 3D segmentation was first delineated automatically by means of a fixed threshold of 40% of the SUVmax, which were corrected by a radiologist manually afterward, blinded to surgical and pathological results. We adopted three steps to preprocess the PET and CT images prior to feature extraction (17). Firstly, we resampled all images to a uniform voxel size of 1 mm × 1 mm × 1 mm using linear interpolation to minimize the influence of different layer thicknesses. Secondly, based on the gray-scale discretization process (bin width for CT = 25, bin width for PET = 0.1), we convert the continuous image into discrete values. Finally, we use the Laplacian of Gaussian and wavelet image filters to eliminate the mixed noise in the image digitization process in order to obtain low-or high-frequency features. Radiomics features were extracted from each PET-derived volume of interest (VOI) and CT-derived VOI by applying dedicated AK software (Artificial Intelligence Kit; GE Healthcare), which is in compliance with image biomarker standardization initiative guidelines (18). A total of 2074 radiomics features were extracted from each VOIs (1037 for CT, 1037 for PET) including (i) 198 for first-order feature, (ii) 14 for shape feature, (iii) 264 for gray level co-occurrence matrix (GLCM) feature, (iv) 176 for gray level size zone matrix (GLSZM) feature, (v) 176 for graGy level run length matrix (GLRLM) feature, (vi) 55 for neighborhood gray tone difference matrix (NGTDM) feature, (vii) 154 for gray level dependence matrix (GLDM) feature.

Feature Selection
After the radiomics features extraction, Z-score normalization was done on each radiomics feature. In addition, the same preprocessing procedure was also applied to the testing set. The dataset was randomly assigned to either the training set or test set in 7:3 ratios. Intra-and inter-class correlation coefficients (ICCs) were calculated to assess the intra-and inter-observer reproducibility, and those radiomics signatures with ICC lower than 0.80 were excluded due to the poor reproducibility. Specifically, Reader 1 and Reader 2 drew the VOIs of 60 cases (20 EGFR-WT NSCLCs, 20 EGFR-19-MT NSCLCs and 20 EGFR-21-MT NSCLCs) of CT images and PET images randomly selected from the whole cohort. Reader 1 repeated the segmentations two weeks later. ICC greater than 0.80 indicated good agreement of feature extraction. The VOI segmentation for the remaining cases were performed by Reader 1.
The feature selection was carried out by using a stepwise selection method. Firstly, univariate logistic regression analysis was utilized to select features with P < 0.05 for the subsequent analysis. Secondly, multivariate logistic regression analysis was applied to choose features closely related to different EGFR status. The P-in and P-out of multivariate logistic analysis were 0.05 and 0.10, respectively. Finally, a subset of the most  informative features was retained using the least absolute shrinkage and selection operator (LASSO) method.

Machine Learning Model
Based on clinical variables, PET metabolic parameters, and PET/ CT-derived radiomics features, three different machine learning classifiers were applied to develop a comprehensive model for differentiating between EGFR-WT, EGFR-19-MT, and EGFR-21-MT, respectively. A support vector machine (SVM) model was built bused on the selected optimal feature subsets of the training dataset. The hyper-parameters of the SVM model were automatically selected by the search method. The kernel, gamma and C were "rbf", 0.1 and 0.1, respectively. Similarly, two other models using RF and DT classifiers were also established.

Construction of Radiomics Nomograms
For patients receiving TKIs targeted therapy, all the clinical prognostic factors (including EGFR mutation site, gender, smoking status, pathological stage, location, histologic subtype, CEA, age and tumor size) and PET metabolic parameters (SUVmax, SUVmean and TLG) were evaluated by univariate analysis using the Kaplan-Meier approach. Statistically significant variables were analyzed for the multivariate Cox forward stepwise regression model to select independent predictors of OS and PFS. Cox regression models were utilized to select the most useful predictive features associated with patients' survival outcomes. A PET/CT radiomics score (Radscores) was calculated for each patient by a linear combination of selected features weighted according to their respective coefficients, and corresponding nomograms were established by integrating the independent prognostic indicators as well as the Rad-score to assess survival benefit. To assess the clinical usefulness of the nomograms, C-index was calculated to evaluate the performance of the models, calibration curve analysis and DCA were performed for estimating and comparing the clinical usefulness of nomograms.

Treatment, Follow Up and Survival Analysis
All patients with EGFR mutation type received first-line EGFR-TKI therapy and routine follow-up after treatment. The endpoints of this study were PFS and OS. PFS is defined as the time interval from treatment to recurrence or progression of the disease. OS is defined as the time interval from treatment to death. Survival curves were drawn using the Kaplan-Meier approach and compared using the log-rank test. Censored data were removed and all remaining data were used for survival analysis.

Statistical Analysis
Univariate analysis (chi-square test or Mann-Whitney U test) was performed by using SPSS software (Version 25.0, IBM). The predictive performance of the machine learning models was determined by the receiver operating characteristic (ROC) curve, and area under the curve (AUC) were calculated. The "RMS" package was used to create the nomogram (19). All statistical analyses of this study were performed using R 3.5.1 and Python 3.5.6. A double-tailed P value less than 0.001 indicated statistical significance.

Intra and Inter-Observer Reproducibility of Feature Extraction
The intra-observer ICC ranged from 0.809 to 0.914, and interobserver ICC ranged from 0.758 to 0.900, therefore, an ideal intra-and inter-observer reproducibility of feature extraction was demonstrated in our study.

Feature Extraction and Selection
A total of 2632 radiomics features were extracted from each VOIs (1316 for CT, 1316 for PET), and 14 radiomics features were filtered, which consisted of six CT-derived radiomics features and eight PET-derived radiomics features. The radiomic features and corresponding coefficients are listed in Supplementary Table 3.

Performance of Different Prediction Models
The ROC analysis demonstrated clinical usefulness of the SVM model, which is superior to the DT model and RF model. All results regarding diagnostic efficacy were displayed in Table 2 and the ROC curves were demonstrated in Figure 3.

Construction and Validation of Radiomics Nomogram
Among clinical parameters, SUVmax and mutation sites proved to be independent predictors of OS and PFS, which was integrated into the nomogram's development in Supplementary Tables 4-7.
Radiomics features for calculating PET/CT Rad-scores of OS and their importance were displayed in Table 3. Radiomics features for calculating PET/CT Rad-scores of PFS and their importance were displayed in Table 4.  The C-index of the integrated nomogram in the training and validation sets were 0.80 and 0.82, respectively. The integrated nomogram outperformed the radiomics nomogram and the clinical nomogram. Nomograms were shown in Figure 4. The diagnostic performance of nomograms is shown in Table 5.
The corresponding calibration curve and decision curve are displayed in Figures 5, 6.

DISCUSSION
In summary, there are two highlights of our study. Firstly, we developed the first-of-its-kind PET/CT-derived radiomic signature based on the three-classification approach, which demonstrated excellent clinical usefulness in predicting EGFR mutation status. The radiomic signature successfully stratified   In this study, we firstly explored the potential association between PET metabolic parameters and the EGFR mutation profiles. Our findings demonstrated that there was a significant difference in SUVmax between EGFR-WT, -19-MT and -21-MT patients. Similarly, in a previous study conducted by Lv et al. (20) confirmed that 18 F-FDG PET/CT metabolic parameters' values were significantly lower in EGFR-MT than in EGFR-WT NSCLCs. Another previous study also reported that EGFR-MT lung adenocarcinomas have relatively lower 18 F-FDG uptake in comparison with EGFR-WT tumors (21), and SUVmax of patients EGFR-21-MT was higher than that of EGFR-19-MT (22). The possible reasons are explained as follows: EGFR mutation was correlated with low tumor metabolic activity of NSCLCs on 18 F-FDG PET/CT. Several researchers considered that EGFR-TKIs could accelerate the glucose uptake of tumor cells. Specifically, tumor cells with high glucose metabolism levels have abundant glucose uptake. Thus, they have less demand for EGFR-TKIs compared to low metabolic tumor cells. As a result, the incidence of EGFR-MT in NSCLCs with high SUVmax is relatively lower (23). Our results are in line with such conclusions. However, different from other acceptable notions, Results from Lee et al. (24) and Minamimoto et al. (25). Indicated that no significant difference was found regarding the SUVmax between the EGFR-WT and EGFR-MT patients, suggesting that SUVmax was not an independent predictor for EGFR mutation. Previous studies conducted by Kanmaza et al. (26) and Ko et al. (27) demonstrated that a higher SUVmax was associated with an EGFR mutation. As a result, these conflicting results demonstrated that 18 F-FDG uptakes may not be a dependable marker for predicting EGFR mutation status. The possible reasons for these discrepant findings can be attributed by the patient baseline demographics of the enrolled patients, the small study sample size number of patients in our study, and the complex tumor microenvironment.
Although a significant relationship between the tumor glucose metabolism level depicted on PET images and EGFR mutation profiles has been reported in several works (22,28,29), traditional PET-derived semiquantitative indexes show insufficient ability to be widely used in clinical practice. It has been demonstrated that SUVmax as a single pixel value only yield moderate AUC for differentiating EGFR-WT from EGFR-MT, whereas TLG as a volumetric measurement of glucose metabolism level has not demonstrated more satisfactory performance either. Thus, our study established a comprehensive prediction model based on 18 F-FDG PET/CT radiomics analysis to provide additional value in optimizing the predictive performance for EGFR mutation profiles in patients with NSCLC.
Radiomics, as an emerging field, has greatly promoted the diagnostic and prognostic accuracy. Currently, radiomics for determining gene mutation status in patients with NSCLC based on PET/CT images were reported in several studies (29)(30)(31). In a previous study, Zhang et al. (32) developed radiomics model to assess the predictive power of pre-therapy 18 F-FDG PET/CTbased radiomic features for EGFR mutation status in NSCLC. However, firstly, it was carried out on a relatively small sample size (two hundred and forty-eight patients). Secondly, the area under the curve values analysis for predicting EGFR mutation status displayed limited discrimination performances (with AUC equal to 0.79 in the training set, and 0.85 in the validation set). Thirdly, advanced radiomics features were not extracted for all patients for technical reasons (only 47 PET and 45 CT radiomic features). In contrast, multiple machine learning classifiers were utilized to identify predictive radiomic features, and the SVM model yielded a training AUC of 0.881, 0.851 and 0.849 in EGFR-WT, EGFR-19-WT and EGFR-21-WT, respectively, whereas a validation AUC of 0.926, 0.805 and 0.859, respectively in the current study, which might provide higher diagnostic performance.
The current study was applied relatively larger sample size, higher-order features and advanced radiomics analysis methods, as well as high-dimensional radiomics signatures extracted up to 2632. Li et al. (33) developed radiomics model through an integrated analysis of 115 NSCLC patients with somatic mutation testing to investigate the feasibility of quantitative and qualitative features extracted from PET-CT in evaluating EGFR mutation status in NSCLC patients. Only a total of 38 radiomic features quantifying tumor morphological, grayscale statistic, and texture features were extracted from the primary PET/CT images. A radiomic signature based on both PET and CT radiomic features outperformed individual radiomic features, the PET or CT radiomic signature. Additionally, a combined radiomic signature with clinical factors exhibited a further improved performance in EGFR mutation status differentiation in NSCLC. In the present study, we also constructed the different classifiers based on integrated radiomic features derived from PET, CT and metabolic parameters to further improve the diagnostic ability. Furthermore, in terms of predicting EGFR gene mutations in NSCLC, few studies involve predicting the certain EGFR mutation site (EGFR-19-MT or EGFR-21-MT) using PET-CT. The study of Zhang et al. (34) have validated that only one PET radiomics feature demonstrated significant but low predictive ability (AUC = 0.661) for differentiating EGFR-19-MT from EGFR-21-MT. Compared with the above study, our prediction model can distinguish EGFR-WT, EGFR-19-MT and EGFR-21-MT in one stop, and shows good discrimination performance.  Regarding strengths of the present work, our results not only predicted EGFR mutation status and mutation site, but also predicted patient survival outcomes, which have scarcely been investigated. In clinical practice, although the tumor, node, and metastasis (TNM) staging system are commonly applied to evaluate the survival prognosis of malignant tumors, we have to admit that this method still has many inevitable shortcomings in the prognostic assessment of lung cancer (35). In fact, the survival period of patients at the same stage may differ. Thus, a one-size-fits-all strategy based on TNM is not applicable in all situations. Novel methods of prognostic assessment are urgently needed to achieve precision treatment. In this study, we supplied clinicians with an easy-to-use method for predicting survival outcomes in NSCLC patients receiving targeted therapy by constructing a radiomics nomogram that exhibited excellent performance, with high c-indexes in the validation set. Furthermore, with the inclusion of clinic-pathological variables in a single nomogram, the prediction performance was further improved, which could allow for better decision-making for NSCLC patients. In the present study, we found that SUVmax and mutation site were independent predictors of the survival period, suggesting their clinical usefulness in the long-term management of NSCLC patients receiving TKIs. Our data provided concordant results to previous study that SUVmax can provide some evidence for survival prognosis (36). On the basis of this fact, we guess that the higher the level of glucose metabolism, the more aggressive tumor cell growth is, and the poorer the patient's survival prognosis is (37). Yang et al. demonstrated that gender was an important prognostic risk factor in NSCLC patients receiving TKI therapy, which is inconsistent with our findings, possibly due to differences in the inclusion of the study population (38). Although this study has obtained satisfactory results, there are still several limitations: Firstly, patient selection bias might exist due to the retrospective nature. Thus, a prospective validation might provide sufficient evidence for further clinical application. Secondly, cases from a single center and relatively small sample size may impair the portability of the prediction model. It is necessary to conduct multi-center research to enhance the generalization ability of the model. Thirdly, only lung adenocarcinoma was included in this study. The predictive ability of our model in other lung cancer types is needed to be validated. Fourth, as for the delineation of the lesions, a semiautomatic segmentation method is used. The more timeconsuming approach should be explored in the future.

CONCLUSION
In conclusion, this study demonstrated that the pre-treatment PET/CT-based radiomics features exhibited excellent performance for the prediction of EGFR mutation profiles in lung adenocarcinoma. Furthermore, we provided an easy-to-use approach to predict the survival outcome of patients receiving targeted therapy, which can be very useful in the clinical practice to guide individualized molecular targeted therapy.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Harbin Medical University Cancer Hospital. The ethics committee waived the requirement of written informed consent for participation.


Adaptor linked K63 di-ubiquitin activates Nedd4/Rsp5 E3 ligase

Nedd4/Rsp5 family E3 ligases mediate numerous cellular processes, many of which require the E3 ligase to interact with PY motif containing adaptor proteins. Several arrestin-related trafficking adaptors (ARTs) of Rsp5 were self-ubiquitinated for activation, but the regulation mechanism remains elusive. Remarkably, we demonstrate that Art1, Art4, and Art5 undergo K63-linked di-ubiquitination by Rsp5. This modification enhances the plasma membrane recruitment of Rsp5 by Art1 or Art5 upon substrate induction, required for cargo protein ubiquitination. In agreement with these observations, we find that di-ubiquitin strengthens the interaction between the pombe orthologs of Rsp5 and Art1, Pub1, and Any1. Furthermore, we discover that the homologous to E6AP C-terminus (HECT) domain exosite protects the K63-linked di-ubiquitin on the adaptors from cleavage by the deubiquitination enzyme Ubp2. Together, our study uncovers a novel ubiquitination modification implemented by Rsp5 adaptor proteins, underscoring the regulatory mechanism of how adaptor proteins control the recruitment, and activity of Rsp5 for the turnover of membrane proteins.

Introduction 24
The Nedd4/Rsp5 family E3 ligases are responsible for membrane protein ubiquitination, required 25 for endocytosis and lysosome-dependent protein degradation. Tryptophan-tryptophan (WW) 26 domains of Nedd4 family E3 ligases bind to substrate proteins via interaction with PY motifs 27 containing a consensus sequence P/L-P-x-Y (Rotin & Kumar, 2009;Schild et al, 1996). Other 28 substrates lack PY motifs and instead rely on interactions with adaptor proteins that recruit the 29 Nedd4 E3 ligase to them, exemplified by a family of arrestin-related trafficking adaptors (ARTs) 30 that bridge the association between substrates and Rsp5 for ubiquitination (Lin et al, 2008). Since ARTs ubiquitination enhances function, we next sought to test if C-terminal fusion with 184 ubiquitin molecules could rescue the cargo sorting defect of art1 K486R or art1 ∆PY . Since toxic 185 arginine analog canavanine is transported by PM transporter Can1 in yeast and Can1 endocytosis 186 prevents subsequent cell death (Grenson et al, 1966), canavanine hypersensitivity occurs when 187 Can1 cannot be endocytosed (such as in an art1∆ mutant), which provides a readout of Art1 188 function. Thus, we examined the canavanine sensitivity of the art1∆ mutant expressing Art1 WT , 189 art1 K486R or art1 ∆PY fused with C-terminal 1x, 2x or 3xUb. The C-terminal fusions with ubiquitin 190 molecules did not enhance the functionality of Art1 WT , art1 K486R or art1 ∆PY (Fig. S2I, S2J). Besides 191 Art1, we tested if Itr1-GFP sorting can be restored by art5 K364R with C-terminal 1xUb or 2xUb and 192 found the 1x or 2xUb fusions do not enhance the Itr1 sorting (Fig. 2F). Together, our data indicate 193 that di-ubiquitin needs to be conjugated at specific residues for proper functionality. The art5 K364R mutant partially blocks the ubiquitination and cargo sorting of Itr1 after inositol 197 treatment (Figure 2A and 2D), but still interacts with Itr1. We therefore hypothesize that the 198 defective ubiquitination of art5 K364R may impair Rsp5 recruitment to the PM. To test this idea, we 199 examined the localization of Art5-GFP in yeast cells before and after inositol treatment. As seen 200 in figure 3A, the Art5 WT -GFP localized at cytosol, nucleus and occasional cytosolic puncta (Sec7-201 negative, Figure 3A). Strikingly, the Art5 WT -GFP is re-localized to PM puncta and patch structures 202 after 1 hour of inositol (20µg/ml) treatment. This result is in line with our previous finding that 203 Art5 is partially translocated to peripheral puncta after shift from minimal media to YPD (Lin et 204 al., 2008), probably because the inositol concentration in YPD is not high enough to drive 205 significant Art5 re-localization and Itr1 endocytosis. In comparison to Art5 WT , the art5 K364R -GFP 206 and art5 ∆PY mainly remain in the cytosol even after inositol treatment ( Figure 3B-D). We conclude 207 that ubiquitination of Art5 is important for protein re-localization to the PM upon inositol treatment. 208 We next asked if Rsp5 can be re-localized to the PM in an Art5-dependent manner after adding 209 inositol to the growth media. As expected, Rsp5 was observed to be recruited to PM patches after 210 inositol treatment in WT cells. However, Rsp5 PM recruitment after inositol treatment was 211 substantially reduced in cells expressing either art5 K364R or art5 ∆PY (Figure 3E and 3F (Lin et al., 2008). As seen in the figure S3A, Art1 is efficiently recruited 218 to the PM in YPD or in methionine media. In contrast to Art1 WT , the recruitment of art1 K486R to 219 the PM is largely attenuated and no PM recruitment is seen with art1 ∆PY (Figure S3B, S3C) . We 220 next tested whether Art1 facilitates PM recruitment of Rsp5. As expected, methionine treatment 221 induces Rsp5 PM recruitment in cells expressing WT Art1, but this recruitment is much reduced 222 in cells expressing art1 K486R ( Figure S3E, S3F). Taken together, our results support the model that 223 specific ubiquitination of adaptor proteins is required for proper recruitment of Rsp5 to target 224 membranes and subsequent ubiquitin-mediated endocytosis of cargo proteins. 225 226 Substrate dependent PM recruitment of adaptor protein Art5 and Art1 227 Next, we sought to examine if cargo proteins are required for adaptor protein recruitment to their 228 functional locations. To do so, we examined the Art5-GFP PM recruitment in ITR1-WT and itr1∆ 229 mutants upon inositol treatment for 1 hour. Strikingly, we found that the PM recruitment of Art5-230 GFP is abolished in the itr1∆ mutant ( Figure 3G, 3H). Similarly, we observed that the PM 231 recruitment of Art1 is attenuated in the mup1∆ mutant with methionine induction for 1 hour (figure 232 S3G and S3H). We further tested whether Art1 can be recruited to the PM in cells expressing 233 mup1-Q49R mutant upon methionine treatment. Previous data showed that Mup1 mutant Q49R 234 is unable to be endocytosed with methionine treatment (Guiney et al., 2016) and the mup1-Q49R 235 mutation abolishes the protein-protein interaction between Mup1 and Art1 ( Figure S2E). We 236 expected that PM recruitment of Art1 will be impaired in this mutant due to the loss of the 237 interaction between Art1 and Mup1-Q49R. Indeed, PM recruitment of Art1 is abrogated in the 238 mup1-Q49R condition ( Figure S3G and S3H), suggesting that the Mup1-Art1 interaction is 239 required for methionine-induced Art1 PM re-localization. Collectively, our data demonstrate that 240 the substrate proteins (Itr1 and Mup1) are required for adaptor protein (Art5 and Art1) recruitment 241 to target membranes (Fig. 3I). 242

K63-linked di-ubiquitination enhances the interaction between adaptor proteins and Rsp5 244
We found that the Rsp5 adaptor proteins (Art5, Art1, and Art4) undergo K63-linked di-245 ubiquitination and this modification is required for efficient recruitment of Rsp5 to target 246 membranes and cargo protein degradation. We hypothesized that adaptor di-ubiquitination 247 enhances protein-protein interactions between di-ubiquitinated adaptors and Rsp5 and thus 248 promotes the recruitment of the E3 ligase. To test this hypothesis, we set out to examine the binding 249 between mono-Ub or K63-linked di-Ub and the HECT domain of Rsp5, as well as the binding 250 between PY motifs and WW domains. To do so, we first generated K63-linked ubiquitin chains showed that M1 linked di-Ub binds to HECT with Kd=36µM, comparable with the K63-linked di-283 Ub (Fig. 4D). In line with this result, our in-vitro binding assay showed that the binding between 284 GST-2xUb and HECT domain is stronger than GST-Ub (Fig. 4E). In comparison, K48-linked di-285 Ub shows a much lower affinity than K63-di-Ub, Kd=145µM (Fig. 4F). Together, our results 286 demonstrate that HECT domain specifically binds to linear form K63 di-Ub and the exosite site is 287 required for ubiquitin binding. di-Ub cannot be made (Tsui et al, 2005). We found that the K63-linked Ub I44A -Ub WT binds to 296 HECT with an Kd=120µM, lower binding affinity than the K63 di-Ub (Fig. 4G). Thus, our result 297 suggests that both distal and proximal ubiquitins contribute to the HECT domain binding, probably 298 cooperatively. 299

300
We next sought to determine if K63-linked di-Ub enhances the binding between adaptor and 301 HECT type E3 ligase. In spite of the fact that KR mutants of Art5 and Art1 lead to attenuated 302 Rsp5 PM recruitment and cargo proteins (Itr1 and Mup1), we still observed the interaction between 303 Art1-K486R or Art5-K364R with Rsp5 using Co-IP (Fig. 4H, 4I), probably due to the interaction 304 between the PY motifs and WW domains. Indeed, as shown in the figure (Fig. 4J, 4K), Art1 and 305 Art5 PY motif containing peptides interact with purified WW domains from Rsp5 (Kd= 3.6µM 306 for Art1 PY motifs and Kd=3.1µM for Art5 PY motifs), but not with PY motif mutants. Since the 307 sole interaction between PY motifs and WW domains does not suffice the full activation of Rsp5 308 function ( Fig. 2A and S2A), the interaction between di-Ub and HECT may enhance the binding 309 affinity between Need4/Rsp5 E3 ligases and their adaptors. We next sought to test the binding 310 between full length adaptors and Rsp5. We found that we could not express Art1 or Art5 at high 311 levels in E. coli, then tried to express the Art1 orthologue from S. pombe, Any1, in E.coli. We 312 found the S. pombe Rsp5 ortholog Pub1 interacts with Any1 with a binding affinity Kd~2.1µM 313  We then sought to uncover the mechanism by which the Art5 degradation is triggered. We 364 observed that Art5 protein is also degraded in the Ub-K63R mutant (Fig. 1C). We wondered if 365 ubp2∆ can rescue Art5 degradation in the Ub-K63R mutant. To test this, we deleted Ubp2 in the 366 Ub-K63R mutant and found that ubp2∆ does not reverse the loss of Art5 protein (Fig. 5D). Our 367 results suggest that the loss of K63-di-Ub on Art5, instead of Ubp2, in either the Ub-K63R mutant 368 or rsp5-exosite mutant, leads to Art5 degradation. Remarkably, in these two conditions, the Art5 369 degradation can be rescued in art5-K364R and art5-∆PY mutant ( Fig. 5B and 5D). These data 370 suggest that Rsp5 can mediate both K63 and K48-linked ubiquitination and the K364 residue of 371 Art5 can be conjugated with both K63-linked di-Ub and K48-linked polyubiquitin chain. 372 Collectively, our results support a working model that the K63-linked di-Ub on Art5 is fully 373 engaged into the Rsp5 exosite so that Ubp2 cannot cleave it efficiently, whereas the K63-linked 374 di-Ub is disengaged in an rsp5-exosite mutant therefore exposed to Ubp2 for cleavage. Upon the 375 cleavage of K63-linked di-Ub, a K48-linked polyubiquitin chain is conjugated at the same residue 376 Art5-K364 thereby leading to the proteasome-dependent degradation of Art5 (Fig. 5E-5H). 377

378
Since PM recruitment of Rsp5 is enhanced by Art5 or Art1 protein ubiquitination (Fig. 3E, S3E) 379 and Ubp2 can deubiquitinate these adaptor proteins ( Fig. 5A and S4D), we therefore asked if the 380 adaptor ubiquitination process is reversible and Ubp2 is involved in this process or not. To monitor 381 the pre-existing adaptor proteins, we decided to employ the tet-Off system to fix the pool of adaptor 382 proteins by treating the cells with doxycycline. As seen in the figure S4F and S4G, pre-existing 383 Art5 or Art1 undergoes ubiquitination upon inositol or methionine treatment for 1 hour in both 384 WT and ubp2∆ conditions, whereas the adaptor proteins shifted back to less ubiquitinated status 385 after removing the stimulation in the WT condition, but not in the ubp2∆ mutant. Together, our 386 data demonstrate a model for adaptor protein recycling mediated by Ubp2 (Fig. S4H). First, 387 stimulation enhances adaptor ubiquitination. Second, the ubiquitinated pool of adaptor proteins 388 can be de-ubiquitinated by Ubp2 when stimulation is terminated.

K63-linked di-Ub is engaged into Rsp5 E3 ligase for activation 409
While we showed that Rsp5 adaptors Art1, Art4 and Art5 undergo K63-linked di-Ub modification, 410 we also demonstrate that this conjugation sharply enhances the binding with the E3 ligase and 411 activates the E3 ligase function for substrate ubiquitination (Fig. 2D). We reason that the 412 interaction between the di-Ub chain and the HECT domain locks the E3 ligase and adaptor into an 413 active/functional conformation. For adaptor-independent ubiquitination, the Nedd4/Rsp5 ligase 414 exosite is also required for efficient ubiquitin conjugation, demonstrating that the "Ub-exosite 415 binding" is required to localize and orient the distal end ubiquitin chain to promote conjugation 416 (Kim et al., 2011;Maspero et al., 2011). In terms of the Rsp5 adaptor-mediated function, we 417 propose that the binding between "di-Ub and exosite" not only enhances the binding affinity 418 between the E3 ligase and adaptor (Fig. 4L-4M), but also leads to more productive Rsp5 419 recruitment to properly orient and present the substrate for ubiquitination at target membranes (Fig.  420   3B). 421

422
While we presented the evidence of E3 ligase activation by ubiquitinated adaptors, we also showed 423 that K63 di-Ub generates a 6-fold tighter binding to the HECT domain than mono-Ub. We reason 424 that the K63 di-Ub provides alternative options to bind a single site, but also fits with a model in of Rsp5 via various adaptors leads us to hypothesize that activated Rsp5 can be released from 477 engaged adaptor proteins. We showed that the adaptor proteins Art1 and Art5 undergo di-478 Ubiquitination upon environmental stimulation and Ubp2 is required to reverse this ubiquitination. 479 Once the ubiquitination is done, the engaged K63 di-Ub is exposed for cleavage by Ubp2. 480 Thereafter, Ubp2 acts on ubiquitinated adaptor proteins to release the adaptor proteins and Rsp5. 481 The mechanism by which Ubp2 executes this reaction needs to be addressed in the future.