Bedside functional monitoring of the dynamic brain connectivity in human neonates

Patient recruitment

Nine healthy very preterm neonates (five males and four females) born at 28 ± 2 weeks PMA were included in the study. The data from the three first patients were used to set up the fUS headset design and therefore were not included in the quantitative study. Four healthy term-born neonates (two males and two females) in their first week of life were also included. Finally, one term newborn with a KCNQ2 mutation was imaged as well. This patient exhibited a pathological “burst-suppression” EEG pattern, associated with frequent ictal seizures. In general, mutations in KCNQ2 encoding the voltage-gated potassium channel Kv7.2 are found not only in 60–70% of families with the autosomal dominant self-limiting syndrome of benign familial neonatal seizures, but are also associated with a severe neonatal epileptic encephalopathy with early-onset epileptic encephalopathy within the first week of life, as reported here³¹. EEG at onset shows a “burst-suppression” pattern⁵ in a majority of these severely affected patients. After the neonatal period, various patterns of seizures have been reported. Cognitive outcome varies widely from mild to moderate intellectual disability to severe developmental delay.

All the acquisitions were approved by the local ethical committee (Robert Debré Hospital, Comité de Protection des Personnes #120601, BELUGA protocol, promoted by INSERM, Institut National de la Santé et de la Recherche Médicale, French Health Institute). They strictly complied with the principles for medical research involving human subjects of the World Medical Association Declaration of Helsinki. Written and informed consent was obtained from all parents guardians.

Ultrasound scanner and neonate headset

The fUS acquisitions were performed at the patient’s bedside by regular ultrasound operators. An ultrafast ultrasound research scanner (Aixplorer, SuperSonic Imagine, Aix-en-Provence, France) was used in conjunction with a custom linear ultrasonic probe (128 elements, 0.2 mm pitch, 6.4 MHz, Vermon, Tours, France) placed on the neonate’s anterior fontanel. As fUS rely on the detection of subtle CBV variations typically during several minutes, it is therefore critical to ensure the stability of the probe on the fontanel. This constraint was addressed by conceiving a custom headset (Fig. 1a, e–g) combining a 3D-printed probe holder and a bio-compatible silicon casting adapted to newborns skull shape (Elite double 8, Zhermack, Badina Polesine, Italy). The headset was gently secured to the baby’s head using a harness of straps. The silicon chamber was filled with ultrasonic gel and the probe could be positioned in a coronal or sagittal configuration. Then, it could be tilted manually or electronically using a small servomotor plugin (TowerPro MG91). The total setup weighed 50 g for the manual version, and 65 g for the motorized version. Installing the headset on the patient typically took <1 min. While imaging in conventional B-Mode, the position of the cap was manually adjusted on the fontanel to find the best acoustic window and to ensure the orientation of the acoustic plane. The headset was compatible with non-invasive ventilation system commonly used on very preterm neonates. As a 60 min clinical EEG recording was performed at the same time, the headset remained on the head during this period of time, even though fUS was only acquired during 20 min. This way, fUS acquisition did not interfere with the standard of care. No side-effect associated to fUS acquisition has been reported. The same headset was used for preterm and term neonates. However, as the morphologies of these two populations can be strongly different, the headset was designed to be adaptable to different skull shapes. Thus, for the same probe-holder, different caps can be designed to fit any kind of skull. This setup was used to maintain the probe in a fixed plan during fUS acquisition phase but also to acquire 3D ultrasound volume during the tomography phase (see Fig. 1f, g for details).

EEG hardware and windowing

In order to make sure that fUS observations arise from comparable neurological conditions, the wakefulness or sleep state of the patients as well as some physiological parameters (cardiac rhythm, respiration, and movements) were monitored. We used a Nihon-Kohden EEG-1200 system to record the EEG, electrocardiogram (ECG), and respiration. We used a standard eight-electrodes longitudinal bipolar montage in the 10/20 convention, placed around the ultrasound headset. The EEG device was synchronized to the ultrasound scanner using a TTL protocol. The preterm and healthy term newborns were observed during phases of quiet sleep. This sleep state was easily identifiable thanks to well-known EEG patterns such as the “trace discontinue”⁵. Quiet sleep was chosen because of patients’ calmness, which limits the occurrence of motions artifacts. Typically, these phases last 20 min. As for the neonate with abnormal development, we focused on the burst-suppression phases due to their long duration (10–20 min) and to their clinical interest.

Power doppler images

Ultrasound can provide vascular information through the so-called Doppler modes. Particularly, the power doppler of a given pixel of the image has been known to be proportional to the fractional moving blood volume in this pixel. In the context of ultrafast ultrasound, the power doppler was acquired as follows. The tilted ultrasonic plane waves (−3°, 0°, and 3°), were successively emitted by the probe with a pulse repetition frequency of 1.8 kHz (PRF). Their backscattered echoes were coherently added and reconstructed into an image through a post-processing operation called beamforming. This resulting compounded image was therefore obtained at a frame rate of 600 Hz. A stack of 342 of such compounded images was acquired in 570 ms, followed by a pause of 430 ms enabling data processing and probe cooling. For this 1-second long elementary block, the tissue clutter was suppressed using an adaptive singular value decomposition (SVD) based filter^32,33. This automatic SVD process ensure optimal and unbiased blood signal filtering. The remaining blood signal was integrated over the whole 342 frame stack, giving rise to a high-definition power doppler image. Such image was acquired in every scanned plane for the tomography phase or in a given plane every 1 s during 20 min for the fUS imaging phase. This sequence was far below the FDA requirements in terms of acoustic power, mechanical and thermal indexes (ISPTA = 30 mW/cm², MI = 0.5, TI = 0.5).

The SNR of the power doppler images is directly related to the number of plane-waves used for coherent compounding. Should the SNR be too low, it would be possible to increase the number of tilted plane waves and the PRF in each elementary block. It has to be noted that the PRF and the imaging depth were capped by hardware constraints. The greater the number of frames and the imaging depth, the higher the data amount to be processed for each block. The data pipeline has to be optimized to allow high PRF and deep images. This optimization work was still in progress when the patient with burst-suppression showed up. This is why for this patient, the imaging depth was limited to 35 mm (vs 49 mm for other patients, after optimization). Recent scanners and computers have now overcome this limitation.

Ultrasound tomography

Using the optional motorized stage of the probe holder, a plane-by-plane scan could be automatically conducted in coronal or sagittal views. In each plane, both the ultrafast power doppler and the conventional anatomic B-Mode were acquired (Fig. S1). Knowing the 3D cylindrical coordinates of each plane, the plane-by-plane scan was interpolated on 3D cubic grid to form a volume. Once the motorized sagittal scan was over, the probe had to be manually rotated of 90 degree to perform the coronal scan. This operation was facilitated by the magnets (Q-05-1.5-01-N, Supermagnete, Gottmadingen, Germany) between the cap and the probe-holder. The cap remained on the patient’s head, whereas the probe-holder was easily detached, rotated, and re-attached. As this manual rotation could cause slight cap displacements, the Doppler volumes obtained with the coronal and sagittal scans were automatically registered using a Matlab routine (imregister). The resulting geometric rigid transformation (T_{Tomo}^{Rigid}) was applied to the B-Mode volumes as well. Doppler volume tubular structures were then enhanced with a Hessian-based filter³⁴ for better 3D rendering, namely, (V_{Doppler}^{Tomo}).

Furthermore, the B-Mode images exhibited well visible hypoechoic structures filled with CSF, namely, the subarachnoid space, the lateral cerebral ventricles and the cavum septum pellucidum when it existed (Fig. S2). These structures can be easily segmented in any coronal and sagittal plane by many different image-processing techniques. We used the combination of two morphological reconstructions in Matlab (imimposemin and imreconstruct). This 2D segmentation allowed the 3D reconstruction of the CSF-filled structures (V_{CSF}^{US}), which was used as a base for the MRI atlas registration step.

MRI atlas registration and region extraction

We use the “Consistent high-definition spatiotemporal neonatal brain atlas” from Makropoulos et al.¹³, which was freely available from the brain Development project website (https://brain-development.org/brain-atlases/neonatal-brain-atlases/neonatal-brain-atlas-makropoulos/). This MRI-based neonatal atlas was chosen among several others²² for three reasons: first and foremost, it was one of the few atlases available for very preterm infants born at 28 weeks PMA; second, it was established on a rather high number of participant at this age (N = 5); third, it came with a great number of parcels (N = 87). The CSF-filled structures were labeled in this atlas and gave access to the volume(V_{CSF}^{MRI}). The two volumes to be registered (V_{CSF}^{US}) and (V_{CSF}^{MRI}) were first coarsely manually aligned to minimize the convergence time of the subsequent registration program. Then, we used a Matlab imregister routine to finely align the volumes and obtain the full geometric transformation (T_{US to MRI}). In order to compensate for small differences in shapes between the MRI atlas and the patient’s anatomy, (T_{US to MRI}) was set as non-rigid (affine transformation comprising translation, scale, shear, and rotation). This allows a stretching of the different structures to improve the matching of the two volumes. This transformation could be also applied to the Doppler volume (V_{Doppler}^{Tomo}) (Fig. 1c) in order to visualize the accessible cerebral volume. Finally, for any ultrasound plane, the transformation (T_{US to MRI}) gave access to the cerebral structures intersected by the ultrasound beam. Thus, these regions’ boundaries were extracted and overlaid on the ultrasound plane to serve as a parcellation (Movie S1). Owing to inter-patient variability, the parcellation could be slightly manually adjusted to avoid aberrant segmentation: typically, this was to ensure that venous sinuses and cortical areas were not mixed up.

It should be noted that the registration performed here is applied to structures segmented on 3D ultrasound volumes, and to these same structures segmented in a MRI atlas. The inputs of the registration program are then binary volumes (voxels of one inside the structure of interest, and 0 outside), and not raw MRI or ultrasound data.

Functional ultrasound signal pre-processing

The fUS sequence is based on the repetition of the 1-second long power doppler elementary block described previously. As power doppler is proportional to the CBV, the fUS sequence gives then access to CBV variation along time. For each 1s-long block, only one cardiac cycle was retained to compute the power doppler image, in order to limit the potential influence of cardiac pace variation. After 20 min, when 1200 frames were obtained, a 2D rigid spatial registration was performed to align every frame with the first one. This was done to take into account the potential slight drifting of the cap. Besides, occasional patient brisk movements could cause cap motion on the skin and result in artifacts. The better the headset is secured to the head, the less frequently these artifacts occur. Therefore, it is technically possible to strongly reduce the occurrence of motions artifact simply by ensuring a tight headset setup. The remaining motion artifacts were automatically identified in the signal and the corresponding time points were rejected (Fig. S6). The CBV signal was then spatially smoothed with a Gaussian kernel whose full-width at half-maximum (FWHM) was 0.5 mm. The resulting signals were found to have dominant spontaneous fluctuations in the < 0.1Hz frequency band. This domain has been extensively described to be associated with resting-state FC dynamics¹¹. Therefore, the CBV was filtered in the [0.01 0.1]Hz frequency band using a bandpass fifth-order Butterworth zero-phase filter.

Static connectivity: correlation matrices

For each cerebral parcel in the ultrasound plane, the pre-processed CBV was averaged over all pixels. For the area #i, we note this temporal signal CBV_i(t), with (overline {{mathrm{CBV}}_i}) being its temporal average and σ_i its standard deviation. For a given number of parcels N_parc and a given patient, the correlation matrix was computed over all the pairs of Pearson correlation coefficients as follow:

Homotopic connectivity: seed-based correlation and homotopic transform

Seed-based correlation maps were computed by manually selecting a pixel and correlating its pre-processed CBV signal with all other pixels signals on 5-min-long windows. The resulting correlation map was thresholded by applying the Fisher’s z-transform and choosing a level of significance of z > 3.72 (P < 0.0001, one-tailed test), corresponding to r > 0.40. The homotopic transformation ({T}_{{R} to {L}}^{Homo}), which links a pixel to its contralateral counterpart was obtained by automatically registering the mask of the right hemisphere Msk_R on the mask of the left hemisphere Msk_L. These masks were obtained by simply fusing the parcels of the FL and the cingulate gyrus. As the midline was not perfectly centered, Msk_R and Msk_L were cropped to have a common surface. Then, ({T}_{{R} to {L}}) was fitted in Matlab by combining a reflection and an affine geometric transform. The PMHC was therefore obtained by computing the correlation for each pixel ({mathrm{Pix}}) as:

Where corr stands for the Pearson correlation. The resolution of ({T}_{{R} to {L}}) was measured by considering the mismatch between ({Msk}_{L}) and ({T}_{{R} to {L}}({Msk}_{R})). This error was typically 0.5 mm, hence Gaussian spatial smoothing FWHM of the CBV was enlarged to 1 mm. Finally, in order to average the PMHC maps between patients, the different morphologies were aligned to a common reference. The first patient was arbitrarily chosen as the reference. The parcellation map of each patient was registered to the reference’s parcellation using imregister routine with affine transforms. These transforms were applied to the corresponding PMHC maps, which were then averaged.

Dynamic connectivity: phase matrices and clustering

As correlation coefficients were only significant for sufficiently long period, we chose another estimator for assessing dFC: the instantaneous phase shift between the different regions’ CBV time series. The implicit assumption was that the real signal CBV_i(t) could be written as ({mathrm{CBV}}_ileft( t right) = a_ileft( t right){mathrm{cos}}left( {varphi _ileft( t right)} right)). Therefore, it is common to access the phase (varphi _ileft( t right)) using the Hilbert transform (H_ileft( t right)) of CBV_i(t):

with j being the imaginary unit. At each time point, the phase shift cosine between two regions’ time series gave a measure of their mutual synchronicity. Typically, a coefficient of one meant that the regions were in-phase, and 0 that they were out-of phase. For two given parcel p and q, we could then define the coefficient m_t(p, q) of the so-called phase matrix M_t at time point t by:

Consequently, a fUS acquisition of N_t time points and N_parc parcels bound to a stack of phase matrices of size [N_parc, N_parc, N_t].

A full stack of phase matrices was obtained over N_acq independent fUS acquisitions by concatenating all the M_t matrices in the third dimension, leading to a stack M^all of size [N_parc, N_parc, N_acq × N_t]. Here, we kept only N_parc = 4 parcels, with N_t = 300 seconds and N_acq = 11 windows (five acquisitions on preterm, three on term, and three on the pathological neonate). It boiled down to a stack of 3300 phases matrices.

Supervised clustering was performed on the complete M^all stack using the Matlab implementations of the k-means algorithm (500 trials, random initialization, L₁-norm distance) with a parametric number of states N_states. The L₁–norm has been reported to be the most robust for this task (17, 18) and was kept for coherence with the literature. No significant changes in interpretation were observed when using L₂-norm though. For each state, the significance of the occurrence rates differences between the three populations was assessed by a one-way ANOVA followed by pairwise post hoc Tukey’s range tests.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.