Referee Consensus: A new enabling technology for computed tomography

1. Referee Consensus: A new enabling technology for computed tomography Don Krieger : Clinical Neurophysiology Malcolm McNeil: Language Disorders Jinyin Zhang, Xin Li: Electrical Engineering Walter Schneider: Psychology David O. Okonkwo: Neurosurgery Catalin Trenchea: Mathematics Bedda Rosario-Rivera: Statistics Jim Becker: Neuropsychology

2. XSede13 #86

3. Support XSEDE: Extreme Computing Consortium OSG: Open Science Grid Army Research Laboratory Pittsburgh Foundation Clinical and Translational Science Institute (CTSI)

4. Our clinical problem is to understand and treat concussion. • We use magnetoencephalography (MEG) during task performance to look for disordered brain activity. • We extract neuroelectric information from the MEG recordings using a novel nonlinear optimization method, referee consensus. • The OSG grid supercomputer was essential to our effort. Summary

5. RESULTS TO DATE:● The method can resolve simultaneously active sources 1-2 mm apart from single trial data

6. RESULTS TO DATE:● Retinal activation is seen.

7. RESULTS TO DATE:● The source count shows task relevant differential activation in 3 mm voxels.

8. Concussion alters thinking: processing speed, cognition, memory, attention. • Neuropsychological measures are sensitive to behavioral correlates of these changes. • MEG measures are sensitive to neuroelectric correlates of these changes. • Clinical objectives: • Identify sites and mechanisms of disordered neurophysiologic function. • Promote development of targeted therapies. Clinical Neuroscience

9. Continuously record the extracranial magnetic field while the volunteer performs a cognitive task … The Measurement

10. A detectable MEG source: electrical current flow within localized neural populations (105 – 108 cells). • The field strength measured at each sensor is a linear combination of fields due to … • an unknown number of sources • at unknown locations (3 values) • with unknown vector amplitude (2 values) • The field strength at the sensors is linear in the amplitudes and nonlinear in the xyz coordinates. The Problem:Neuroelectric source identification

11. ● The measurements are a multivariate time series, i.e.: multiple observations in time ... … at multiple sensors ...● The measurement sensors could be … … photon counters (imaging) … accelerometers (seismic tomography) … sound transducers (ultrasound imaging, SONAR) … direct brain recording electrodes (connectivity analysis) … virtual recording sites (functional connectome) ● Each source must contribute to the observations at many sensors. The Solution:Referee consensus nonlinear optimization

12. ● Calculation of the referee template is a spatio temporal curve fitting procedure and the referee metric quantifies the “goodness of fit.”● The referee consensus cost function is computed in the search space … “measurement error” is not used.● It is insensitive to the presence of sources which are not included in the equations, i.e. it can be used to decide if there is a source present at a specific location regardless of the presence of other sources. The Solution: Referee consensus nonlinear optimization

13. ● Short data epochs (15-100 msec) from single trial data. • ● Comparable to the “gold standard” ECD Localization • 1. Choose a location – this linearizes the model. • 2. Solve the model for the remain linear variables. • 3. Apply the solution to the sequence of MEG observations. • 4. Compute the referee consensus cost function. • 5. If the decision is that a source is present, record it. • 6. If not, use the gradient of the cost function to choose the next location and continue at step 2. The Solution:Referee consensus nonlinear optimization

14. m1 m2 . . . m306 Forward Solution For each point in time, there are 306 MEG measurements.

15. h1,1s1 = m1 m2 . . . m306 Forward Solution The simplest model is that there is a single source with amplitudes1. The forward solution model, h1,1, is a nonlinear function of the location of the source and of the location, orientation, and shape of the sensor at which m1 is obtained. h1,1 is the first entry in the lead field matrix.

16. h1,1s1 = m1 m2 . . . m306 Forward Solution Notice thath1,1s1 = m1is solvableonly whenh1,1is a number, i.e. when the location of the source and the location, orientation, and shape of the sensor atwhichm1 is obtained are fixed.

17. h1,1s1 = m1 m2 . . . m306 Forward Solution Instead of magnetic field detectors, the sensors could be photon counters (imaging), accelerometers (seismology), microphones (ultrasound or SONAR), ... The functional forms of the forward solutions would be different. And so the h ‘s would be different. But the form and substance of these linearized equations and their solution would be the same, i.e. referee consensus is applicable.

18. h1,1s1 = m1 + ε1 m2 . . . m306 Forward Solution We add an error term, ε1, to account for error in accounting for the measurement with a single source and the fixed values used to compute h1,1.

19. h1,1s1 = m1 + ε1 h2,1s1 = m2+ ε2 . . . m306 Forward Solution We add terms for the 2nd equation. The source amplitude, s1, is the same. But the lead field matrix entry, h2,1, is different since the location, orientation, and shape of the sensor at which m2 is obtained is different than for m1.

20. h1,1s1 = m1+ ε1 h2,1s1 = m2+ ε2 . . . . . . h306,1s1 = m306+ ε306 Forward Solution We continue to add terms for the rest of the equations. This formulation is equivalent to that used for Equivalent Current Dipole source localization. For that, the source location is sought which minimizes Σ(ε2), the sum of the squared errors.

21. h1,1s1 = m1+ ε1 h2,1s1 = m2+ ε2 . . h30,1s1 = m30+ ε30 Forward Solution - skip The success of Equivalent Current Dipole source localization relies on the assumption that the measurements are due to a single source. To insure this, the strength of many of the sources is attenuated:1. Averaging to a time locked event is used.2. The measurements are limited to those from the sensors nearest to the source.

22. h1,1s1 + h1,2s2 = m1+ ε1 h2,1s1 + h2,2s2 = m2+ ε2 . . . . . . h306,1s1 + h306,2s2 = m306+ ε306 Forward Solution Continuing, a 2nd source may be added. The 2 source model is the minimum model required for referee consensus. One of the sources is at the location being tested. The other is a dummy location, a referee.

23. h1,1s1 + h1,2s2 + . . . + h1,isi + . . . + h1,nsn + . . . = m1+ ε1 h2,1s1 + h2,2s2+ . . . + h2,isi + . . . + h2,nsn + . . . = m2+ ε2 . . . . . . h306,1s1 + h306,2s2+ . . . + h306,isi+ . . . + h306,nsn + . . . = m306+ ε306 Forward Solution … and many more may be added to the model. If the number of sources exceeds the number of measurements, the solution is “under determined.” Finding a unique solution then becomes dependent on a clever choice of error metric as used in MNE, LORETA, and VESTAL.

24. h1,1s1 + h1,2s2 + . . . + h1,isi + . . . + h1,nsn + . . . = m1+ ε1 h2,1s1 + h2,2s2+ . . . + h2,isi + . . . + h2,nsn + . . . = m2 + ε2 . . . . . . h306,1s1 + h306,2s2+ . . . + h306,isi+ . . . + h306,nsn + . . . = m306+ ε306 Forward Solution If the number of sources is much less than the number of measurements, the solution is “over determined.” The unique solution may be found using ordinary least squares regression. Referee consensus requires this.

25. h1,1s1 + h1,2s2 + . . . + h1,isi + . . . + h1,nsn + . . . = m1+ ε1 h2,1s1 + h2,2s2+ . . . + h2,isi + . . . + h2,nsn + . . . = m2+ ε2 . . . . . . h306,1s1 + h306,2s2+ . . . + h306,isi+ . . . + h306,nsn + . . . = m306+ ε306 Inverse Solution In matrix form the ordinary least squares solution is:S = (HTH)-1HTM S is the set of source amplitudes, the s‘s. H is the lead field matrix. M is the set of measurements, the m‘s.

26. In matrix form the ordinary least squares solution is:S = (HTH)-1HTM Inverse Solution The equation means: “each source amplitude is a weighted sum of the measurements.” The weights are rows of (HTH)-1HT. These rows are called “spatial filters.”1. Each row of (HTH)-1HT has gain of 1.0 for contributions to the MEG from the source corresponding to that row.2. Each row has gain of 0.0 for contributions to the MEG from sources corresponding to the other rows. A formal proof is included in the manuscript. This property is what is used to compute the referee consensus metric.

27. How It Works:The referee template, the activity at the test location, A, from the viewpoint of a referee. There are 4 sources present here including one at A.An additional dummy location is selected, R a “referee.” In this case, there is no source present at R.

28. How It Works:The referee template, the activity at the test location, A, from the viewpoint of a referee. 1. We compute a model due to two sources, A and R. The filter for R applied to the sequence of measurements produces the waveform shown in blue, R!A. It contains nothing of the waveform at A since the filter at R has gain 0.0 for that location. 2. We compute a 2nd model due to two sources, A’ and R. This filter for R produces the similar waveform shown in red, R!A’. But this one contains something of the waveform at location A.

29. How It Works:The referee metric: “Is there a source present at A from the viewpoint of referee R?” The difference between the two filtered waveforms, R!A- R!A’, the lower amplitude waveform in red, is the estimate of the current at location A from the viewpoint of referee R. This is the referee template. This waveform is approximately proportional to the original waveform at location A shown below it.

30. How It Works:The referee metric: “Is there a source present at A from the viewpoint of referee R?” We now compute the quantity of the difference waveform present in each of the two filtered waveforms, i.e. ((R!A- R!A’)● R!A’)2 vs ((R!A- R!A’)● R!A)2. If there is more found in R!A’, the referee metric is set to 1, i.e. from the viewpoint of referee R, there is a source present at A. Otherwise the metric is set to 0.

31. ● 90 referee locations were used. Since the current amplitude has 2 components, there are 2 rows for each location in (HTH)-1HT and 2 referee metrics.● The referee template was generated for each and a consensus template was produced from all 90 x 2 templates. ● The 90 x 2 referee metrics were summed to generate the referee consensus. How It Works

32. ● In order to (1) test for the presence of a source at A and (2) compute the gradient of the referee consensus to decide on the next location to test, the referee consensus was computed 6 times, once each for location A’±1 mm distant along the x,y,z axes. ● 114 of 180 (p < 10-2) was used as the cut-off for each of the 6 tests producing a joint p < 10-12 for acceptance. How It Works

33. Referee consensus is a novel cost function:●The measurements must be linear in at least one variable for each source. Once a location is fixed, i.e. the system is linearized. The “zero” property of ordinary least squares regression is then used to produce estimates of the waveform at the test location from the “viewpoints” of many dummy locations, “referees.”● A sequence of measurements in time must be used. A single snapshot is not sufficient.● Calculation of the referee template is therefore a spatio temporal curve fitting procedure and the referee metric quantified the “goodness of fit.” Use of referee consensus to test for the presence of a source at a specific location: Summary Points

34. Referee consensus is a novel cost function:● Referee consensus is computed in the search space rather than in the measurement space.● It is insensitive to the presence of sources which are not included in the equations, i.e. it can be used to decide if there is a source present at a specific location regardless of the presence of other sources.● This independence enables very efficient implementation on a grid computing resource. Use of referee consensus to test for the presence of a source at a specific location: Summary Points

35. ● Goodness of fit improves with relaxed low-pass filtering and with removal of line noise from the raw MEG.● Simultaneously occurring sources as near as 1 mm to each other are uncorrelated, i.e. they are resolved.● Retinal activation is seen during a visual task.● The source amplitude shows task relevant differential activation in 8 mm voxels.● The source count shows task relevant differential activation in 3 mm voxels.● The task decision impacts subsequent brain activity more strongly in concussed vs normal volunteers. Results to date

36. Goodness of Fit - skip The distribution of referee consensus values: ● improves with less restrictive low-pass filtering, ● This is opposite to what would be expected with a method which uses an error metric ● It suggests that there is useful information present in the MEG above the low pass cut-off, in this case 150 Hz. ● improves with removal of line noise from the raw MEG.

37. Sensitivity to Line Noise - skip ●The averaged fourier transform of the virtual recordings demonstrates contamination by line noise (blue). ● Removal of 60 Hz and harmonics, 120, 180, 240, 300 Hz from the raw MEG with 1 Hz wide notches eliminates this (red). ● Note that the spectrum amplitude (red) from 60-300 Hz is ≈ 1/f, i.e. the power spectrum ≈ 1/f1/2.

38. Resolving power: 1-2 mm Correlation between pairs of 80 msec waveforms for sources 1 to 100 mm apart. Single trial data. This is better than what is achievable with electrocorticography, i.e. direct population recordings.

39. Visualization of Retinal Activation These sources are from a single subject from 120 single trial time segments. Each spanned 40 msec beginning 30 msec after visual stimulation. Note the apparently tonic activation of the cerebellum. 3D Animation (YouTube) 3D Animation (local)

40. Localization of task relevant activation ● Volunteers performed a choice task with continuous MEG recording …referee consensus processing epochs are shown. • Activate many areas of the brain: • Language • Graphic • Decision. • Many trials • Many timing marks

41. The

42. green

43. circle

44. is

45. below.

46. The

47. The

48. The

49. Two trials with event markers … - Normal time scale: about 4 seconds/trial - Stimuli (words and test figures): black arrows - Responses (button presses): magenta arrows

50. Differential brain activation: functional mapping The count of sources found in each 3 mm voxel following a “match” vs a “mismatch” were compared using a X2 statistic. Results are shown only for voxels with p < 10-7.