Wavelet brain angiography suggests arteriovenous pulse wave phase locking

When a stroke volume of arterial blood arrives to the brain, the total blood volume in the bony cranium must remain constant as the proportions of arterial and venous blood vary, and by the end of the cardiac cycle an equivalent volume of venous blood must have been ejected. I hypothesize the brain to support this process by an extraluminally mediated exchange of information between its arterial and venous circulations. To test this I introduce wavelet angiography methods to resolve single moving vascular pulse waves (PWs) in the brain while simultaneously measuring brain pulse motion. The wavelet methods require angiographic data acquired at significantly faster rate than cardiac frequency. I obtained these data in humans from brain surface optical angiograms at craniotomy and in piglets from ultrasound angiograms via cranial window. I exploit angiographic time of flight to resolve arterial from venous circulation. Initial wavelet reconstruction proved unsatisfactory because of angiographic motion alias from brain pulse motion. Testing with numerically simulated cerebral angiograms enabled the development of a vascular PW cine imaging method based on cross-correlated wavelets of mixed high frequency and high temporal resolution respectively to attenuate frequency and motion alias. Applied to the human and piglet data, the method resolves individual arterial and venous PWs and finds them to be phase locked each with separate phase relations to brain pulse motion. This is consistent with arterial and venous PW coordination mediated by pulse motion and points to a testable hypothesis of a function of cerebrospinal fluid in the ventricles of the brain.

The videos are then imported in to the computational analysis environment   (Fig B). 28 Visual inspection the visible and NIR angiogram video files recored by the operating 29 microscope reveals that they are not in synchrony when downloaded from the 30 microscope, but instead are temporally misaligned by up to 11 frames (the image 31 frequency is 30 Hz). Their timings must be adjusted to enable the analysis of motion 32 and flow. At the end of the indocyanine green (ICG) run, the surgeon halts the 33 recording by pressing a button on the microscope handle. Given that the microscope 34 magnification is at 5x, the button press causes a notable simultaneous displacement in 35 both the visible and NIR videos that can be used to synchronize them. I wrote a 36 graphical user interface program to manipulate the temporal alignment of the visible 37 and NIR video footage by reference to the sharp field of view displacement cause by the 38 button press, and estimate that the videos are aligned thereby to within about 1 frame 39 Perflutren (Definity, Lantheus Medical Imaging, North Billerica, Massachusetts) is an 48 ultrasound contrast agent that mimics red blood cell rheology, is hemodynamically inert 49 at diagnostic ultrasound frequency and intensity ranges [3,4], and enables quantitative 50 perfusion imaging [5][6][7][8][9][10][11][12]. Evolving methods ultrasound imaging methods with increased 51 sensitivity and speed appear likely to permit a fidelity and sampling rates beyond those 52 employed here [13,14]. 53 Imaging of the phantom and the piglets was at 8.3 MHz using a linear array probe 54 (MySonos 201, Medison) with ultrasound gain dials (overall, near-field, and far-field) 55 were set to 32 on the device. The analog images produced by the device were digitized 56 (DT3162 Variable-Scan Monochrome Frame Grabber, Data Translation, Marlboro, 57 Massachusetts, for early experiments, and Dazzle DVBridge, Pinnacle Systems, 58 Mountain View, CA, for the latter experiments; I verified in vitro that both systems 59 provided equivalent digitized values for the same phantom) and stored as cine image 60 sequences in QuickTime format [2]. The cine images were then imported into the 61 Mathematica environment for analysis by custom programs. concentration-to-signal relationship, as microspheres tend to rise to the surface in a 66 stagnant system, preventing uniform mixing. Therefore, to simulate active mixing more 67 typical of a physiologic system, I built closed-loop and open-loop phantoms. In vitro 68 phantoms for assessing quantitative ultrasound have been employed by others [15].

69
The closed-loop phantom is designed to simulate steady-state flow. Polyvinyl tubing 70 of volume 30 mL was attached to a peristaltic pump programmed to 100 mL per minute, 71 providing a circular flow rate similar to that observed in mammalian arterial circulation. 72 Images were obtained as the perflutren microsphere concentration was systematically   suggesting that tracer kinetic methods are applicable to the ultrasound signals obtained 82 from perflutren bolus passage in the piglet model described here [16].

84
Prior to computational analysis, for each cine ultrasound sequence a polygon was drawn 85 on a representative ultrasound image to separate the intracranial region from the 86 remainder of the anatomy, as illustrated in Fig F. Based on this boundary, bony skull 87 base structures and extracranial vessels were excluded from all analysis. To capture a 88 motion signal, the ultrasound video was inspected, and a site of motion independent of 89 the passage of contrast was identified. In two animals, this site was at the inferior border 90 of the cranial cavity, and in one it was located peripherally in a fissure. In one animal, a 91 site with consistent motion could not be identified, so it is excluded from further 92 analysis because its vascular pulsations could hence be referenced to pulse motion. A 93 polygon was drawn on each motion site, and a motion signal was measured using the 94 same algorithm and software code as applied to the intraoperative human optical data. 95 The coronal images, motion region of interests, and angiographic time intensity curves  Summary data for the cranial window ultrasound angiograms are presented in Table 98 B;  The same method and the same computer source code is employed to measure cardiac 103 frequency (CF) pulse brain motion for the human brain surface optical and piglet 104 cranial window ultrasound data. An angiogram video is inspected for an object or 105 region with particularly clear pulse motion. Its spatial boundaries are marked by a 106 polygon drawn in a graphical user interface on a still image from the video sequence.

107
The pixels within the polygon define the region of interest (ROI) for motion tracking.

108
To reduce motion noise, the video frames are averaged with each frame replaced by 109 the mean of the 5 nearest frames. The spatial displacement is calculated between each 110 such frame pair 6 frames apart. This is selected instead of adjacent frames because too 111 little motion generally occurred between adjacent frames to be consistently detectable. 112 The policies of local frame averaging and of measuring motion across a frame interval of 113 5 were found by trial and error generally to yield satisfactory signal to noise 114 ratios (SNRs) of the motion signal. The method of approximating CF SNR is described 115 below in Section D.

116
For each such frame pair, the ROI in the second frame is displaced between -5 to +5 117 pixels both horizontally and vertically and each such displaced version is paired with 118 the same first frame, giving 11x11=121 total combinations of ROI pairs. The pixel-wise 119 correlation between the ROI of the first frame and each of the 121 displaced ROIs of 120 the second frame are computed, giving a matrix of 121 correlations. The center of mass 121 of these 121 correlations is taken to represent the net displacement vector of the ROI 122 between the two frames. This displacement vector is then calculated for all successive 123 frame pairs, giving n-5 displacement vectors for n frames, to which is prepended 2 124 copies and appended 3 copies of the null displacement vector { 0 ., 0 . } (Fig I).

125
For convenience, the displacement vector is then reduced in dimensionality from a 126 sequence of two-dimensional displacements to a one dimensional motion signal as 127 follows. The two-dimensional displacements are plotted and treated as a scattergram.

128
The covariance matrix of the scattergram is computed and subjected to 129 eigen-decomposition. Then each individual two-dimensional displacement is projected 130 onto the axis of the major eigenvector of the covariance matrix. This gives a sequence of 131 scalar displacements that termed here the motion signal, M (indexing by uniformly 132 sampled time intervals is implied). The SNRs yielded by these methods are shown for 133 the human data in Table A and for the piglet data in Table B.

135
I employ Fourier methods to look for the presence of CF signal relative to noise before 136 analyzing the time-resolved CF signal by wavelet methods. The Fourier analysis of SNR 137 is based on dividing a narrow band signal of interest (CF in this instance) by a 138 broadband background [17]. However in these data there are other signals in the 139 broadband background that should not be counted as noise. These include the 140 respiratory frequency signal and the very low frequency bolus passage signal. I exclude 141 these from the broadband background in generating the denominator for SNR 142 calculations. Instead the SNR is estimated here as the height of the Fourier spectral 143 peak for CF divided by the mean height of 5 peaks starting about 10 steps away in the 144 higher frequency direction. This appears to serve as a suitable broadband denominator 145 in all signal data analyzed here, including angiographic and pulse motion data for both 146 humans and piglets. The results, presented in Tables A and B, should be taken as mere 147 approximations but they appear to support the presence in all of the data of a CF 148 signal appropriate for further analysis.

151
For a complex-valued CF datum c = a + ib in a spatiotemporal grid, its magnitude |c| is 152 rendered in a pixel as brightness and its phase φ c as hue (Fig H).

153
The complex continuous wavelet transform of f (t) is given by the equation where s is wavelet scale, as above, u is the wavelet translation parameter, and the 155 superscript * represents complex conjugation [18]. The mother wavelet, ψ, represents

176
In this paper, transformation by a high temporal resolution wavelet ψ, symbolized 177 by the over hat symbolˆ, is reified by the library function GaborW avelet [1], and 178 transformation by a high frequency resolution wavelet ψ, symbolized by the over tilde 179 symbol˜, is reified with equal performance by the library functions GaborW avelet [6] 180 and GaborW avelet [12]. The real and imaginary components of GaborW avelet [1] and 181 GaborW avelet [6] are depicted in Fig K. 182 The source code for the wavelet functions is not released by the Mathematica library, 183 but computational testing of these functions against test data discloses that  Table D. For two one-dimensional time signals g(t) and f (t), the cross-correlation h(t) is given by 212 If g(t) and f (t) have wavelet transforms g(t) and f (t), then gives the wavelet transform, h, of the cross-correlation h(t), where as before the 214 superscript * denotes complex conjugation [19,20].

215
Using the notation of this paper (Table C), the proposed reconstruction by wavelet 216 cross-correlation is thus where C * is the one-dimensional motion signal wavelet transformed with a wavelet ψ.

218
The cross-correlated result x has a hat over symbol rather than the unspecified wavelet 219 resolution as x to specify that at inverse wavelet transformation a high temporal 220 resolution wavelet ψ will be employed (with filtering for s ♥ ).

221
Optionally, C may be replaced by the normalized to enhance comparability between subjects and species.

223
The wavelet cross-correlation of Equation E may be simulated as in the method of  Even though the PW images reconstructed from piglet ultrasound have a more 260 grainy appearance than those from human optical angiography, there is greater cardiac 261 frequency signal in the piglet ultrasound images. This may seen by comparing the 262 cardiac frequency SNRs in the piglet ultrasound versus human optical angiographic data 263 (compare SNRs in Tables A and B). The bases of these differences remain unknown.

G Raw Data And Wavelet Angiography Video Files 265
All raw data and computationally reconstruted video data are supplied as Supporting

266
Information. The video file names for the raw data and for the key vascular PW cine 267 images are presented in Table E.   In the spatiotemporal grid of a reconstruction of x, a cumulative distribution function is 298 constructed of the magnitude of the complex-valued grid elements. The shape of the 299 curves shows that about 20% of the grid elements have significant CF signal for both 300 the human and piglet subjects. For constructing the phase histograms the phases are 301 drawn from the grid elements in the upper 20% of CF magnitude (Fig Q). Separate 302 phase histograms are made of the arterial and venous classified pixels, then the two are 303 merged into one (Fig R).         S7 Video S8 Video S9 Video S10 Video S11 Video S12 Video h3 S13 Video S14 Video S15 Video S16 Video S17 Video S18 Video h4 S19 Video S20 GaborWavelet [6] Real Imaginary

Fig K. Temporal And Frequency Resolution Of A Wavelet ψ Function.
Vertical axis is magnitude and horizontal axis is time. The narrow temporal span of the real and imaginary components of the upper row corresponding to the function call in this wavelet computational library of GaborW avelet[1] reflects a high temporal but low frequency resolution of that wavelet ψ. The broad temporal span corresponding to the function call GaborW avelet [6] reflects a low temporal but high frequency resolution. The ridge at s ♥ is preserved in the high frequency resolution (HFR as computed by GaborW avelet [12]) ψ scalograms, but not in the high temporal resolution ψ scalograms, which indicates the presence of frequency aliasing. This simulation shows that use of a pure HTR wavelet ψ may lead to excessive frequency alias in reconstructed images. In the employed wavelet computational library, the function GaborW avelet [1] is a high temporal resolution wavelet transform and GaborW avelet [12] a high frequency resolution transform. Graph (c) is the cross-correlation of (a) and (b). The magnitude and phase filtered by cardiac wavelet scale are in the second and third rows. This one-dimensional simulation suggests that cross-correlation of a high temporal resolution wavelet transform by a high frequency resolution wavelet transform can attenuate some of the frequency alias that accompanies a high temporal resolution wavelet transform alone. A CDF is made of the magnitudes and only those φ values with | • | greater than the 80 th percentile are included in the phase histograms. This reduces the contribution by phase values in low signal pixels. This figure is for h1. Those for the h2-4 and p1-3 are similar in shape.