# Using Graph Theory to Assess the Interaction between Cerebral Function, Brain Hemodynamics, and Systemic Variables in Premature Infants

Authors: Dries Hendrikx; Liesbeth Thewissen; Anne Smits; Gunnar Naulaers; Karel Allegaert; Sabine Van Huffel; Alexander Caicedo
Journal: Complexity
 (2018)

---

## Abstract

Graphs can be used to describe a great variety of real-world situations and have therefore been used extensively in different fields. In the present analysis, we use graphs to study the interaction between cerebral function, brain hemodynamics, and systemic variables in premature neonates. We used data from a propofol dose-finding and pharmacodynamics study as a model in order to evaluate the performance of the graph measures to monitor signal interactions. Concomitant measurements of heart rate, mean arterial blood pressure, arterial oxygen saturation, regional cerebral oxygen saturation—measured by means of near-infrared spectroscopy—and electroencephalography were performed in 22 neonates undergoing INSURE (intubation, surfactant administration, and extubation). The graphs used to study the interaction between these signal modalities were constructed using the RBF kernel. Results indicate that propofol induces a decrease in the signal interaction up to 90 minutes after propofol administration, which is consistent with clinical observations published previously. The clinical recovery phase is mainly determined by the EEG dynamics, which were observed to recover much slower compared to the other modalities. In addition, we found a more pronounced loss in cerebral-systemic interactions with increasing propofol dose.


---

## Body

## 1. Introduction

A graph is a structure that can be used to represent the relation between different objects. In this context, a graph can be thought of as a diagram which consists of a set of points, where some or all of them are joined by lines. Formally, the points of the graph are referred to as vertices or nodes, whereas the lines between them are called edges or links. In general, graphs can be used to describe a great variety of real-world situations [1]. Think, for example, of a social network, where people are represented by nodes and the edges between the nodes are used to indicate friendship. Another example is a geographic network of cities, with an edge between two cities indicating a direct connection through a highway. In addition to the presence (or lack) of an edge connecting two nodes, extra measurements can be associated with the edges. These measurements are formally referred to as edge weights. In a social network, edge weights could be used to denote the strength of the friendship (acquaintances, close friends, …). In a geographic network, the weights can indicate the physical distance or the amount of traffic typically encountered on each road. Mathematically, this type of diagram corresponds to a weighted graph.In the present analysis, weighted graphs are used to study the interaction of cerebral function, brain hemodynamics, and systemic variables in premature neonates. Multiple studies are available in the literature that studied the pairwise interactions between some of these variables. Caicedo et al. analysed the relation between mean arterial blood pressure (MABP) and regional cerebral oxygen saturation (rScO2), measured by means of near-infrared spectroscopy (NIRS) [2]. The coupling between these two variables, defined using a transfer function approach, was found to be a measure to assess cerebral autoregulation. Semenova et al. examined the relation between MABP and electroencephalography (EEG) [3]. The authors documented that preterm infants with a high clinical risk index for babies (CRIB) score were found to be associated with a higher nonlinear coupling between EEG activity and MABP, quantified by means of mutual information. Tataranno et al. examined the relation between rScO2 and EEG and found that increased oxygen extraction was related to spontaneous activity transients observed in the EEG [4]. In contrast to the studies mentioned above, we aim to analyse the interaction between cerebral and systemic variables using an extended multimodal approach, integrating three systemic variables: heart rate (HR), MABP and arterial oxygen saturation (SaO2), rScO2, and EEG.This study is situated within the interdisciplinary field of network physiology, which analyses how diverse physiologic systems dynamically interact and collectively behave to produce distinct physiologic states and functions [5]. Moreover, the use of graphs enables a graphical representation of the interaction between the different physiological systems in time. This study shows for the first time a comprehensive model of different physiological processes comprising autoregulation, neurovascular coupling, or baroreflex, working at the same moment in time. In literature, most studies focus on these processes individually without taking into account the influence of the other processes. With the graph approach outlined in this paper, we try to show the different processes, their interaction, and the importance of the individual processes at each moment in time. To the best of our knowledge, this is a totally new mindset and way of showing the physiological interaction between cerebral function, brain hemodynamics, and systemic variables in newborn neonates.The interaction between the different variables is studied using premedication by means of propofol as a model. Propofol (2,6 diisopropylphenol) is a short-acting anesthetic: it has a rapid onset of action and is generally short in duration. In neonates, however, it is documented that clinical recovery takes time [6]. In clinical practice, propofol is administered to the neonates as a single intravenous (IV) bolus. Propofol administration is frequently associated with a decrease in MABP in neonates [6–11], children [12], and adults [13–15]. Propofol distributes into the central nervous system and fat tissue immediately after intravenous dosing, which explains the rapid onset of this anesthetic drug. In a secondary phase, propofol is redistributed into the circulation, which leads to vasodilation. Combined with the blunted reflex tachycardia, this can result in hypotension [10]. Therefore, a decrease in MABP is observed up to one hour after administration of propofol in neonates [8]. Premedicating neonates with propofol generally causes a modest and short-lasting decrease in HR, SaO2, and rScO2, as opposed to the longer-lasting and more pronounced decrease in MABP [8], [11, 16, 17]. In addition, the discontinuity pattern of the EEG is also influenced by propofol, which induces a reversible state of diminished responsiveness behaviorally similar to quiet (nonrapid eye movement (NREM)) sleep [18]. During quiet sleep, the EEG of premature neonates shows a spontaneous, physiological discontinuity of electrical activity, characterized by higher amplitude, lower-frequency EEG rhythms (tracé alternant (TA)) [19, 20]. This phenomenon is generally referred to as burst suppression, which corresponds to an increase in interburst interval (IBI) duration [21, 22]. Moreover, a larger IBI duration is associated with smaller FTOE values, which indicate lower brain energy consumption [23].This paper is structured as follows. Section2 describes the dataset used in the present analysis. Section 3 discusses the methods, which include EEG processing, the construction of the graph models, and the definition of features computed from the graph models to quantify the strength of the effect of propofol on these interactions. Section 4 presents the results of the paper, which are extensively discussed in Section 5. Finally, Section 6 summarizes the conclusions.


## 2. Dataset

The dataset used in the present analysis was collected as part of a study on propofol dose selection by Smits et al. [6]. In the study, 50 neonates were sedated using propofol as part of an endotracheal intubation procedure. All subjects in the group of study were recruited at the NICU of the University Hospitals Leuven, Gasthuisberg. The trial was registered on ClinicalTrials.gov NCT01621373, and ethical approval was provided by the ethical committee at the University Hospitals Leuven.Due to incomplete data and overly noisy channels found in 28 neonates, only 22 of the 50 neonates are included in this study. These neonates were all sedated using propofol as part of an INSURE (intubation, surfactant, and extubation) procedure. The neonates are characterized by median (range) postmenstrual ages (PMA) of 30 (26–35) weeks and a median (range) dose of propofol (Diprivan 1%; AstraZeneca, Brussels, Belgium) of 1.0 (0.5–4.5) mg·kg−1. In the present analysis, the neonates are stratified into three groups, based on PMA, since this is a major covariate of propofol clearance in the absence of variability in postnatal age (PNA) [24]. These groups are generally referred to as extremely preterm (group 1: <28 weeks PMA), very preterm (group 2: 28–31 6/7 weeks PMA), and moderate to late preterm (group 3: 32–36 6/7 weeks PMA) [25]. Most of the neonates have a PNA of 1 day. For details regarding the composition of the patient groups, the PNA of the patients, and the doses of propofol administered to the subjects of each group, see Table 1. More information about the clinical characteristics of the subjects can be found in the original paper by Smits et al. [6].Table 1
Stratification of the neonates into three age groups, based on postmenstrual age (PMA) in weeks. For each group, the number of patients, postnatal age (PNA) of the patients, and propofol dose values administered to the subjects in the group are presented.
Group 1
Group 2
Group 3
Extremely preterm
Very preterm
Moderate to late preterm
<28 weeks
28–31 6/7 weeks
32–36 6/7 weeks
Number of patients
5
13
4
PNA (days)
Median
1
1
1
Range
1-2
1-2
1–3
Dose (mg·kg−1)
Median
1.0
1.0
2.75
Range
0.5–1.5
0.5–2.5
1.0–4.5Practices on propofol dosing, particularly in highly vulnerable premature neonates, are not standardized and vary between different NICUs. Multiple studies, however, indicate that propofol dose values of 2.0 to 2.5 mg·kg−1 should be used as preintubation medication in premature neonates [9–11, 16]. The dataset used in the present analysis was collected with the aim to find the median effective dose (ED50) of propofol for sedation. Therefore, lower values of starting propofol dose were used, as indicated in Table 1. More specifically, administered dose ranges from 0.5 to 4.5 mg·kg−1 [6]. In general, the oldest neonates were sedated using higher propofol doses compared to the youngest neonates, as can be observed from Table 1.The multimodal dataset used in this study consists of concomitant measurements of five signal modalities, comprising HR, MABP, SaO2, rScO2, and EEG, recorded from 5 minutes before propofol administration up to 10 hours after. For each neonate, a 6-hour long segment of multimodal data was considered in the analysis, where t = 0 was aligned with the moment of propofol administration. This length was defined based on the shortest recording found in the dataset. Thus, all signals were shortened to six hours for all patients in order to provide uniformity. Moreover, the use of a long time window of 6 hours allows focusing on the regime of interest, since we can study the effect of propofol together with the recovery of the neonates from the drug. Propofol is a three-compartment drug, characterized by a short α and β (median estimates of 1 and 13 minutes, resp.) and a long γ half-life (median estimate of 350 minutes) [26, 27]. The pharmacodynamic effects are primarily associated with the first (α) and second (β) exponential half-life, which indicates that the effect of propofol at the end of the analysis window is minimal. This is confirmed by Smits et al., Vanderhaegen et al., and Ghanta et al., who all observed a clinical recuperation from single intravenous bolus propofol administration within the first hour in neonates [6, 8, 16]. Therefore, the analysis window is divided in two parts: the first 3-hour long time window is used to study the response of the neonates to propofol and the intubation procedure, while the last 3 hours are used as reference. Figure 1 presents an example of a 6-hour long segment of multimodal data for one neonate from the group of study.Figure 1
Illustration of the 5 signal modalities used to construct graph models for the neonates. A 6-hour long segment of multimodal data is presented for one neonate in the group of study (PMA 27 weeks, 0.5 mg·kg−1).The systemic variables (HR (beats/min), MABP (mmHg), and SaO2 (%)) were measured with IntelliVue MP70 (Philips, Eindhoven, The Netherlands) with a Nellcor pulse oximeter. These variables were recorded continuously with a sampling frequency of 1 Hz (Rugloop; Demed, Temse, Belgium). All 22 neonates incorporated in the present analysis had an arterial line, which enabled an invasive measurement of MABP. NIRS was used to measure rScO2 (%) noninvasively with INVOS 5100 using a cerebral neonatal OxyAlert NIRS sensor (Covidien, Mansfield, Massachusetts). As for the systemic variables, the sampling frequency for rScO2 is equal to 1 Hz. Cerebral functioning was assessed using a one-channel EEG (μV). The EEG was measured between the C3 and C4 electrodes according to the international 10–20 system with a sampling frequency of 100 Hz (Olympic Cerebral Function Monitor 6000, Natus). EEG segments with impedance values exceeding 10 kΩ were removed from the raw EEG signal [28]. In addition, movement artifacts identified as rapid changes in the impedance measurement were detected and also removed from the raw EEG signal.


## 3. Methods

### 3.1. Running Interburst Interval Duration

In general, EEG signals of premature neonates alternate between periods of activity, called bursts or burst intervals (BIs), and periods of suppressed activity, referred to as IBIs. Thus, the morphology of neonatal EEG is discontinuous, as indicated by the IBIs. However, this discontinuous pattern evolves towards a more continuous trace with increasing PMA. Therefore, some studies have investigated the use of the length of the IBIs as a marker for maturation [29, 30].Due to the different temporal characteristics between the EEG and all other signal modalities, the EEG signals are processed in order to obtain surrogates for brain activity in a similar time frame as the other measured signals. The EEG signal is segmented in burst and IBI segments using an in-house algorithm based on the line length [31]. The root mean squared (RMS) value and the duration in time for burst and IBIs in overlapping windows of two minutes are used as a surrogate for EEG. The running window is shifted in one second, producing a new score every second. In this way, the sampling frequency of the surrogate measures for EEG has the same sampling frequency as the other signal modalities.In total, five features are computed from the discontinuous neonatal EEG: running RMS values of the original EEG, BIs, and IBIs and running duration values of the BIs and IBIs. In this paper, we only report the results using the running IBI duration, since this is a very robust measure for EEG activity, and thus cerebral metabolism, as validated by our group in a previous study [31]. In addition, this measure is highly interpretable. It is important to note, however, that the other EEG features indicate similar results, since the different feature values are highly related. An example of the five EEG features is presented in Figure 2.Figure 2
Illustration of the features computed from the EEG signal. (a) illustrates a 6-hour long EEG segment for one neonate in the group of study (PMA 27 weeks, 0.5 mg·kg−1). (b) illustrates the running RMS value. (c) and (d) illustrate the running RMS and running duration values for BIs (black) and IBIs (gray), respectively.


### 3.2. Graph Model Developed for This Study

In order to quantify the common dynamics of the different signal modalities, and changes thereof due to propofol, the interaction between the variables is modeled using a graph, as illustrated in Figure3. In general, a graph is defined by a nonzero number of vertices (nodes) and a number of edges (links, connections) between these nodes. The model for the neonates is constructed using a complete graph. A complete graph is characterized by the presence of an edge between all the vertices. The vertex set V of the graph consists of n=5 vertices, corresponding to the 5 signal modalities measured in the present analysis, that is,
(1)V=vHR,vMABP,vSaO2,vrScO2,vEEG.Figure 3
Physiological network representing the interaction between 5 signal modalities recorded on a neonate after propofol administration. The graph consists of 5 vertices, corresponding to the signal modalities. In addition, an edge is present between every pair of nodes (complete graph). Each edge is defined by a weight value that represents the interaction between the corresponding signal modalities.A complete graph withn vertices has m=nn−1/2 edges. Therefore, the edge set E of the graph considered here consists of 10 edges. The vertices of the graph model defined in (1) are connected by edges. These edges are defined by the corresponding edge weight values, which are generally used to assess the strength of the connection between a pair of vertices.The topology of the complete graph described in (1) is assumed to be fixed in time. The edge weights, however, change in time, which we hypothesize to reflect the changes in the interaction between the different signals. In order to compute the graph models, the signals are first normalized to N(0,1), since we are interested in the assessment of common dynamics (signal trends in time) and not absolute values of the signals. Next, the edge weights are computed using a 15-minute long running window of multimodal data, which is shifted by 1 minute (14 minutes overlap). Thus, new edge weight values are computed every minute. Finally, two types of interaction curves are extracted from the graph models: the pairwise interaction between two signal modalities, represented by the time course of the corresponding edge weight, and the overall signal interaction, represented by the graph average degree (see Section 3.4).In the present analysis, weight values are used to denote the interaction between two vertices, that is, two signal modalities. If two modalities are characterized by common nonlinear interactions, they follow the same trends in time. We compute the pairwise similarity using two different similarity measures. Consequently, we generate two graph models for each neonate. Both similarity measures use the radial basis function (RBF) kernel, which is a nonlinear similarity measure. As such, the similarity of the different signals is assessed in a possibly infinitely dimensional feature space, defined by the nonlinear mapϕ. However, the similarity in this feature space is computed implicitly using the RBF kernel function. The first similarity measure kTxi,xj uses the raw signals in the RBF kernel and is thus defined as
(2)kTxi,xj=ϕxi⊤ϕxj=exp−∥xi−xj∥22σ2,where xi and xj represent two segments of multimodal data [32] (subscript T indicates that time domain signals are used for the Euclidean distance in the exponent of the RBF kernel). In the present analysis, xi and xj are segments with a length of 15 minutes, as mentioned before. The similarity kTxi,xj is bounded by 0 (absence of common interactions) and 1 (exact common interactions). The signal similarity computed by (2) is a function of the Euclidean distance between input signals. Consequently, it highly depends on signal amplitudes and can be affected by delays between the signals. A graph model computed using the similarity measure kTxi,xj is denoted as GT.The second similarity measure uses the power spectral density (PSD) of the signals in the RBF kernel. Thus, the time input data is transformed to the frequency domain, before computing the RBF kernel function. Mathematically, this similarity measurekFxi,xj is defined as
(3)kFxi,xj=exp−Sxi−Sxj22σ2,where Sxi and Sxj represent the PSD of input signals xi and xj (length of 15 minutes), respectively (subscript F indicates that frequency domain signals are used for the Euclidean distance in the exponent of the RBF kernel). The PSD is computed using Welch’s method using overlapping subwindows of 5 minutes in order to reduce the noise in the PSD estimate (with use of Hamming window, overlap of 4 minutes and 59 seconds). Note that the kernel presented in (3) is a valid positive definite kernel, since the input data is transformed before application of the kernel function. As before, the similarity defined by kFxi,xj is bounded by 0 and 1. The transformation to the frequency domain allows to include time-delayed signal interactions and interactions of opposite sign, in contrast to kTxi,xj which only takes into account instantaneous amplitude interactions. In physiological systems, it is possible that if one signal increases (decreases), another signal decreases (increases) to maintain homeostasis and that this interaction is not instantaneous but delayed. A graph model computed using the similarity measure kFxi,xj is denoted as GF.


### 3.3. Kernel Tuning

In order to compute the similarity measure, the bandwidthσ of the RBF kernel should to be tuned, that is, optimized to avoid kernel overfitting and underfitting. In the present analysis, the similarity measures kTxi,xj and kFxi,xj both depend on this parameter σ. The optimization procedure is the same for both similarity measures. Therefore, it is outlined in terms of kxi,xj, which represents the two similarity measures. The strategy used to select the kernel bandwidth for the present analysis considers kernel matrix Ω, which is defined as
(4)Ωij=kxi,xj,i,j=1,…,n.Note that the kernel matrixΩ is defined by the kernel bandwidth σ through the definitions presented in (2) and (3). The kernel bandwidth σ is tuned by maximizing the Shannon entropy of kernel matrix Ω. The Shannon entropy HΩ is defined as
(5)HΩ=−∑kpklog2pk,where pk is equal to the probability of seeing the kth possible element of matrix Ω. The entropy is thus determined by estimation of the probability density function (PDF) of matrix Ω. By maximizing the Shannon entropy, we try to obtain a uniform distribution of the values in the kernel matrix, and therefore, we avoid overfitting as well as underfitting.The kernel bandwidth is tuned for each neonate individually. The tuned bandwidth is denoted asσopt. The following optimization problem is defined to estimate σopt:
(6)σopt=maxσHΩC,with
(7)ΩC=Ω1Ω2…ΩN,where ΩC is a collection of kernel matrices, computed from all the signal segments recorded per neonate. Thus, a collection of kernel matrices is computed from the 6-hour long data segment instead of only one kernel matrix in the optimization procedure. If we would consider only one kernel matrix per neonate, it would only contain 25 entries, since the kernel matrix is a 5 × 5 matrix. Clearly, this is not enough data to estimate a robust PDF. Therefore, to solve this problem, we assume that the graph model does not change and that it is situated in the same nonlinear subspace throughout the 6-hour long analysis window. This assumption indicates that σopt should be uniform throughout the analysis window and that σopt can be computed using a concatenation of kernel matrices ΩC, as defined in the optimization problem in (6) and (7).Figure4 illustrates the optimization procedure in a schematic way. The original data segment of 6 hours was segmented into nonoverlapping segments of 15 minutes. Thus, N=24 signal segments of 15 minutes were defined. For each of these segments l, kernel matrix Ωl was computed and all these kernel matrices Ωl were concatenated as indicated in (7). The use of a collection of kernel matrices allows to estimate the probability density function, and consequently, the Shannon entropy. Therefore, HΩC is characterized by one global maximum. For the group of study, median (range) values of σopt are 27 (26–29) and 94 (86–113) for kTxi,xj and kFxi,xj, respectively.Figure 4
Method used to tune the kernel bandwidthσ. In (a), the data is segmented in nonoverlapping signal segments of 15 minutes. For all of these segments, a kernel matrix Ωl is computed using a predefined σgs (4). All the individual kernel matrices Ωl are concatenated in ΩC, which is depicted in (b). Next, the Shannon entropy of ΩC is computed. This procedure is repeated for a range of σ values. The σ value associated with maximal HΩC is selected as the bandwidth for the kernel function.


### 3.4. Graph Measures

In order to assess the overall interaction of the multimodal dataset, the average degree of the graph is used. This section introduces the adjacency matrixA of a graph, the degree di of a vertex, and the average degree δG of a graph G.


#### 3.4.1. Adjacency Matrix

A weighted graphG consists of a nonempty finite set V of elements called vertices vi (or nodes) and a finite set E of distinct unordered pairs of distinct elements of V called edges wij (or links) [33]. Note that the edges of the graph are represented by their weights wij. The adjacency matrix A is a matrix commonly used to define the graph G. The adjacency matrix A denotes the presence of edges between the vertices vi of V and their corresponding weights. More precisely, the adjacency matrix A is constructed as
(8)Aij=wij,ifthereisanedgebetweenυiandυj,0,otherwise.


#### 3.4.2. Vertex Degree

The degreedj associated with a vertex vj of an undirected weighted graph G, with adjacency matrix A, is defined as the sum of all edges incident to vj:
(9)dj=∑i=1nwij=∑i=1nAij,where n is the number of vertices. Therefore, the degree dj characterizes the connection strength of the vertex vj with respect to the other vertices of the graph. In practice, the weights of the edges of a graph are often restricted to a predefined range, which is often normalized to wij∈0,1. Considering normalized weights, the degree is bounded by 0 and n−1, where n is the number of vertices of the graph, that is,
(10)0≤dj≤n−1if∀i,j:wij∈0,1.Ifdj = 0, vertex vj is called an isolated vertex, since it is not connected to any other vertex of the graph. A vertex degree dj = n − 1 indicates a dominating vertex vj, connected to all other vertices of the graph with edge weight equal to 1.


#### 3.4.3. Average Degree

The average degreeδG of a graph G is defined as the mean value of all vertex degrees dj
(11)δG=1n∑j=indjand is a measure associated with the overall connectivity of the graph. Evidently, the bounds of δG are equal to those of the individual vertex degree dj defined in (10). Small values (close to 0) imply a weak connectivity, whereas high values (close to n − 1) indicate a very strong connectivity of the graph.


### 3.5. Features to Quantify Interaction Strength

In order to quantify the strength of the changes in signal interaction, two features are computed from the interaction curves: the normalized areaS between the interaction curve and reference level and the maximal deviation ∆ from the reference level. Both feature values are computed in a time frame from 0 to 90 minutes after propofol administration. Reference levels are defined as the median value of an interaction curve in a time frame from 180 to 360 minutes after propofol administration, as mentioned before. Normalization of S is done by dividing the area by the length of the time interval. Note that S and ∆ are bounded by 0 (no deviation from the reference level) and 1 (very strong deviation from the reference level). Figure 5 presents a graphical example of S (Figure 5(a)) and ∆ (Figure 5(b)). The features are computed from the interaction curves in order to assess the effect of propofol on the dynamical interactions among the different signal modalities. In addition, we investigated how these features change with PMA and propofol dose.Figure 5
Features used to quantify the reduction in signal interaction strength:S (gray shaded area) (a) and ∆ (gray arrow) (b). The feature values are illustrated for one neonate in the group of study (PMA 30 weeks, 2.5 mg·kg−1), where the pairwise interaction was computed using kTxi,xj (2). Feature values S and ∆ are computed from 0 to 90 minutes, while the reference level is defined as the median value of the interaction curve from 180 to 360 minutes.
(a)
(b)In the present analysis, the relation between the feature valuesS and ∆ (dependent variables) and PMA and propofol dose (predictor variables) is studied using linear regression models. The coefficient of determination Ri2 is used to indicate the goodness of fit of the linear model (subscript i denotes the predictor variable i). In addition, the coefficient of partial determination was computed to account for the effect of both predictor variables at the same time. The significance of the coefficient of (partial) determination was assessed using the Monte Carlo permutation test with 105 repetitions. A p<0.05 was defined to be statistically significant. A single asterisk, double asterisks, and triple asterisks denote a p value smaller than 0.05, 0.01, and 0.001, respectively.


### 3.6. Implementation

The analysis, the corresponding computations, and figures presented throughout this study are implemented using MATLAB Release 2016b (The MathWorks, Natick, Massachusetts). Graph theory analysis is performed using the MATLAB toolbox for network analysis, provided by MIT Strategic Engineering [34].


## 3.1. Running Interburst Interval Duration

In general, EEG signals of premature neonates alternate between periods of activity, called bursts or burst intervals (BIs), and periods of suppressed activity, referred to as IBIs. Thus, the morphology of neonatal EEG is discontinuous, as indicated by the IBIs. However, this discontinuous pattern evolves towards a more continuous trace with increasing PMA. Therefore, some studies have investigated the use of the length of the IBIs as a marker for maturation [29, 30].Due to the different temporal characteristics between the EEG and all other signal modalities, the EEG signals are processed in order to obtain surrogates for brain activity in a similar time frame as the other measured signals. The EEG signal is segmented in burst and IBI segments using an in-house algorithm based on the line length [31]. The root mean squared (RMS) value and the duration in time for burst and IBIs in overlapping windows of two minutes are used as a surrogate for EEG. The running window is shifted in one second, producing a new score every second. In this way, the sampling frequency of the surrogate measures for EEG has the same sampling frequency as the other signal modalities.In total, five features are computed from the discontinuous neonatal EEG: running RMS values of the original EEG, BIs, and IBIs and running duration values of the BIs and IBIs. In this paper, we only report the results using the running IBI duration, since this is a very robust measure for EEG activity, and thus cerebral metabolism, as validated by our group in a previous study [31]. In addition, this measure is highly interpretable. It is important to note, however, that the other EEG features indicate similar results, since the different feature values are highly related. An example of the five EEG features is presented in Figure 2.Figure 2
Illustration of the features computed from the EEG signal. (a) illustrates a 6-hour long EEG segment for one neonate in the group of study (PMA 27 weeks, 0.5 mg·kg−1). (b) illustrates the running RMS value. (c) and (d) illustrate the running RMS and running duration values for BIs (black) and IBIs (gray), respectively.


## 3.2. Graph Model Developed for This Study

In order to quantify the common dynamics of the different signal modalities, and changes thereof due to propofol, the interaction between the variables is modeled using a graph, as illustrated in Figure3. In general, a graph is defined by a nonzero number of vertices (nodes) and a number of edges (links, connections) between these nodes. The model for the neonates is constructed using a complete graph. A complete graph is characterized by the presence of an edge between all the vertices. The vertex set V of the graph consists of n=5 vertices, corresponding to the 5 signal modalities measured in the present analysis, that is,
(1)V=vHR,vMABP,vSaO2,vrScO2,vEEG.Figure 3
Physiological network representing the interaction between 5 signal modalities recorded on a neonate after propofol administration. The graph consists of 5 vertices, corresponding to the signal modalities. In addition, an edge is present between every pair of nodes (complete graph). Each edge is defined by a weight value that represents the interaction between the corresponding signal modalities.A complete graph withn vertices has m=nn−1/2 edges. Therefore, the edge set E of the graph considered here consists of 10 edges. The vertices of the graph model defined in (1) are connected by edges. These edges are defined by the corresponding edge weight values, which are generally used to assess the strength of the connection between a pair of vertices.The topology of the complete graph described in (1) is assumed to be fixed in time. The edge weights, however, change in time, which we hypothesize to reflect the changes in the interaction between the different signals. In order to compute the graph models, the signals are first normalized to N(0,1), since we are interested in the assessment of common dynamics (signal trends in time) and not absolute values of the signals. Next, the edge weights are computed using a 15-minute long running window of multimodal data, which is shifted by 1 minute (14 minutes overlap). Thus, new edge weight values are computed every minute. Finally, two types of interaction curves are extracted from the graph models: the pairwise interaction between two signal modalities, represented by the time course of the corresponding edge weight, and the overall signal interaction, represented by the graph average degree (see Section 3.4).In the present analysis, weight values are used to denote the interaction between two vertices, that is, two signal modalities. If two modalities are characterized by common nonlinear interactions, they follow the same trends in time. We compute the pairwise similarity using two different similarity measures. Consequently, we generate two graph models for each neonate. Both similarity measures use the radial basis function (RBF) kernel, which is a nonlinear similarity measure. As such, the similarity of the different signals is assessed in a possibly infinitely dimensional feature space, defined by the nonlinear mapϕ. However, the similarity in this feature space is computed implicitly using the RBF kernel function. The first similarity measure kTxi,xj uses the raw signals in the RBF kernel and is thus defined as
(2)kTxi,xj=ϕxi⊤ϕxj=exp−∥xi−xj∥22σ2,where xi and xj represent two segments of multimodal data [32] (subscript T indicates that time domain signals are used for the Euclidean distance in the exponent of the RBF kernel). In the present analysis, xi and xj are segments with a length of 15 minutes, as mentioned before. The similarity kTxi,xj is bounded by 0 (absence of common interactions) and 1 (exact common interactions). The signal similarity computed by (2) is a function of the Euclidean distance between input signals. Consequently, it highly depends on signal amplitudes and can be affected by delays between the signals. A graph model computed using the similarity measure kTxi,xj is denoted as GT.The second similarity measure uses the power spectral density (PSD) of the signals in the RBF kernel. Thus, the time input data is transformed to the frequency domain, before computing the RBF kernel function. Mathematically, this similarity measurekFxi,xj is defined as
(3)kFxi,xj=exp−Sxi−Sxj22σ2,where Sxi and Sxj represent the PSD of input signals xi and xj (length of 15 minutes), respectively (subscript F indicates that frequency domain signals are used for the Euclidean distance in the exponent of the RBF kernel). The PSD is computed using Welch’s method using overlapping subwindows of 5 minutes in order to reduce the noise in the PSD estimate (with use of Hamming window, overlap of 4 minutes and 59 seconds). Note that the kernel presented in (3) is a valid positive definite kernel, since the input data is transformed before application of the kernel function. As before, the similarity defined by kFxi,xj is bounded by 0 and 1. The transformation to the frequency domain allows to include time-delayed signal interactions and interactions of opposite sign, in contrast to kTxi,xj which only takes into account instantaneous amplitude interactions. In physiological systems, it is possible that if one signal increases (decreases), another signal decreases (increases) to maintain homeostasis and that this interaction is not instantaneous but delayed. A graph model computed using the similarity measure kFxi,xj is denoted as GF.


## 3.3. Kernel Tuning

In order to compute the similarity measure, the bandwidthσ of the RBF kernel should to be tuned, that is, optimized to avoid kernel overfitting and underfitting. In the present analysis, the similarity measures kTxi,xj and kFxi,xj both depend on this parameter σ. The optimization procedure is the same for both similarity measures. Therefore, it is outlined in terms of kxi,xj, which represents the two similarity measures. The strategy used to select the kernel bandwidth for the present analysis considers kernel matrix Ω, which is defined as
(4)Ωij=kxi,xj,i,j=1,…,n.Note that the kernel matrixΩ is defined by the kernel bandwidth σ through the definitions presented in (2) and (3). The kernel bandwidth σ is tuned by maximizing the Shannon entropy of kernel matrix Ω. The Shannon entropy HΩ is defined as
(5)HΩ=−∑kpklog2pk,where pk is equal to the probability of seeing the kth possible element of matrix Ω. The entropy is thus determined by estimation of the probability density function (PDF) of matrix Ω. By maximizing the Shannon entropy, we try to obtain a uniform distribution of the values in the kernel matrix, and therefore, we avoid overfitting as well as underfitting.The kernel bandwidth is tuned for each neonate individually. The tuned bandwidth is denoted asσopt. The following optimization problem is defined to estimate σopt:
(6)σopt=maxσHΩC,with
(7)ΩC=Ω1Ω2…ΩN,where ΩC is a collection of kernel matrices, computed from all the signal segments recorded per neonate. Thus, a collection of kernel matrices is computed from the 6-hour long data segment instead of only one kernel matrix in the optimization procedure. If we would consider only one kernel matrix per neonate, it would only contain 25 entries, since the kernel matrix is a 5 × 5 matrix. Clearly, this is not enough data to estimate a robust PDF. Therefore, to solve this problem, we assume that the graph model does not change and that it is situated in the same nonlinear subspace throughout the 6-hour long analysis window. This assumption indicates that σopt should be uniform throughout the analysis window and that σopt can be computed using a concatenation of kernel matrices ΩC, as defined in the optimization problem in (6) and (7).Figure4 illustrates the optimization procedure in a schematic way. The original data segment of 6 hours was segmented into nonoverlapping segments of 15 minutes. Thus, N=24 signal segments of 15 minutes were defined. For each of these segments l, kernel matrix Ωl was computed and all these kernel matrices Ωl were concatenated as indicated in (7). The use of a collection of kernel matrices allows to estimate the probability density function, and consequently, the Shannon entropy. Therefore, HΩC is characterized by one global maximum. For the group of study, median (range) values of σopt are 27 (26–29) and 94 (86–113) for kTxi,xj and kFxi,xj, respectively.Figure 4
Method used to tune the kernel bandwidthσ. In (a), the data is segmented in nonoverlapping signal segments of 15 minutes. For all of these segments, a kernel matrix Ωl is computed using a predefined σgs (4). All the individual kernel matrices Ωl are concatenated in ΩC, which is depicted in (b). Next, the Shannon entropy of ΩC is computed. This procedure is repeated for a range of σ values. The σ value associated with maximal HΩC is selected as the bandwidth for the kernel function.


## 3.4. Graph Measures

In order to assess the overall interaction of the multimodal dataset, the average degree of the graph is used. This section introduces the adjacency matrixA of a graph, the degree di of a vertex, and the average degree δG of a graph G.


### 3.4.1. Adjacency Matrix

A weighted graphG consists of a nonempty finite set V of elements called vertices vi (or nodes) and a finite set E of distinct unordered pairs of distinct elements of V called edges wij (or links) [33]. Note that the edges of the graph are represented by their weights wij. The adjacency matrix A is a matrix commonly used to define the graph G. The adjacency matrix A denotes the presence of edges between the vertices vi of V and their corresponding weights. More precisely, the adjacency matrix A is constructed as
(8)Aij=wij,ifthereisanedgebetweenυiandυj,0,otherwise.


### 3.4.2. Vertex Degree

The degreedj associated with a vertex vj of an undirected weighted graph G, with adjacency matrix A, is defined as the sum of all edges incident to vj:
(9)dj=∑i=1nwij=∑i=1nAij,where n is the number of vertices. Therefore, the degree dj characterizes the connection strength of the vertex vj with respect to the other vertices of the graph. In practice, the weights of the edges of a graph are often restricted to a predefined range, which is often normalized to wij∈0,1. Considering normalized weights, the degree is bounded by 0 and n−1, where n is the number of vertices of the graph, that is,
(10)0≤dj≤n−1if∀i,j:wij∈0,1.Ifdj = 0, vertex vj is called an isolated vertex, since it is not connected to any other vertex of the graph. A vertex degree dj = n − 1 indicates a dominating vertex vj, connected to all other vertices of the graph with edge weight equal to 1.


### 3.4.3. Average Degree

The average degreeδG of a graph G is defined as the mean value of all vertex degrees dj
(11)δG=1n∑j=indjand is a measure associated with the overall connectivity of the graph. Evidently, the bounds of δG are equal to those of the individual vertex degree dj defined in (10). Small values (close to 0) imply a weak connectivity, whereas high values (close to n − 1) indicate a very strong connectivity of the graph.


## 3.4.1. Adjacency Matrix

A weighted graphG consists of a nonempty finite set V of elements called vertices vi (or nodes) and a finite set E of distinct unordered pairs of distinct elements of V called edges wij (or links) [33]. Note that the edges of the graph are represented by their weights wij. The adjacency matrix A is a matrix commonly used to define the graph G. The adjacency matrix A denotes the presence of edges between the vertices vi of V and their corresponding weights. More precisely, the adjacency matrix A is constructed as
(8)Aij=wij,ifthereisanedgebetweenυiandυj,0,otherwise.


## 3.4.2. Vertex Degree

The degreedj associated with a vertex vj of an undirected weighted graph G, with adjacency matrix A, is defined as the sum of all edges incident to vj:
(9)dj=∑i=1nwij=∑i=1nAij,where n is the number of vertices. Therefore, the degree dj characterizes the connection strength of the vertex vj with respect to the other vertices of the graph. In practice, the weights of the edges of a graph are often restricted to a predefined range, which is often normalized to wij∈0,1. Considering normalized weights, the degree is bounded by 0 and n−1, where n is the number of vertices of the graph, that is,
(10)0≤dj≤n−1if∀i,j:wij∈0,1.Ifdj = 0, vertex vj is called an isolated vertex, since it is not connected to any other vertex of the graph. A vertex degree dj = n − 1 indicates a dominating vertex vj, connected to all other vertices of the graph with edge weight equal to 1.


## 3.4.3. Average Degree

The average degreeδG of a graph G is defined as the mean value of all vertex degrees dj
(11)δG=1n∑j=indjand is a measure associated with the overall connectivity of the graph. Evidently, the bounds of δG are equal to those of the individual vertex degree dj defined in (10). Small values (close to 0) imply a weak connectivity, whereas high values (close to n − 1) indicate a very strong connectivity of the graph.


## 3.5. Features to Quantify Interaction Strength

In order to quantify the strength of the changes in signal interaction, two features are computed from the interaction curves: the normalized areaS between the interaction curve and reference level and the maximal deviation ∆ from the reference level. Both feature values are computed in a time frame from 0 to 90 minutes after propofol administration. Reference levels are defined as the median value of an interaction curve in a time frame from 180 to 360 minutes after propofol administration, as mentioned before. Normalization of S is done by dividing the area by the length of the time interval. Note that S and ∆ are bounded by 0 (no deviation from the reference level) and 1 (very strong deviation from the reference level). Figure 5 presents a graphical example of S (Figure 5(a)) and ∆ (Figure 5(b)). The features are computed from the interaction curves in order to assess the effect of propofol on the dynamical interactions among the different signal modalities. In addition, we investigated how these features change with PMA and propofol dose.Figure 5
Features used to quantify the reduction in signal interaction strength:S (gray shaded area) (a) and ∆ (gray arrow) (b). The feature values are illustrated for one neonate in the group of study (PMA 30 weeks, 2.5 mg·kg−1), where the pairwise interaction was computed using kTxi,xj (2). Feature values S and ∆ are computed from 0 to 90 minutes, while the reference level is defined as the median value of the interaction curve from 180 to 360 minutes.
(a)
(b)In the present analysis, the relation between the feature valuesS and ∆ (dependent variables) and PMA and propofol dose (predictor variables) is studied using linear regression models. The coefficient of determination Ri2 is used to indicate the goodness of fit of the linear model (subscript i denotes the predictor variable i). In addition, the coefficient of partial determination was computed to account for the effect of both predictor variables at the same time. The significance of the coefficient of (partial) determination was assessed using the Monte Carlo permutation test with 105 repetitions. A p<0.05 was defined to be statistically significant. A single asterisk, double asterisks, and triple asterisks denote a p value smaller than 0.05, 0.01, and 0.001, respectively.


## 3.6. Implementation

The analysis, the corresponding computations, and figures presented throughout this study are implemented using MATLAB Release 2016b (The MathWorks, Natick, Massachusetts). Graph theory analysis is performed using the MATLAB toolbox for network analysis, provided by MIT Strategic Engineering [34].


## 4. Results

### 4.1. MABP-EEG Pairwise Interaction

The interaction curves of MABP with respect to EEG after administration of propofol att = 0 minutes are illustrated in Figure 6. These curves are computed using kTxi,xj, defined in (2). The EEG signal is represented by the running IBI duration, as outlined in Section 3. From top to bottom, the interaction pattern is shown for the entire group of study (N=22) and the individual age groups presented in Table 1. First, a pronounced loss in interaction is observed, followed by a gradual increase to a reference level, which is in general reached at t = 90 minutes. Note that this loss in interaction is present among all of the signal modalities of the multimodal dataset, as indicated by the graphs in Figure 7. Figure 8 presents the relation between the features used to quantify interaction strength (S and ∆) and PMA and propofol dose. In addition to the data points, the least squares linear fit is defined (straight lines), together with the 0.95 percentiles of the linear fit (shaded area). The goodness of the linear fit is assessed using the coefficient of determination Ri2, which is equal to RA2=0.09 and RD2=0.53 for feature S and RA2=0.17 and RD2=0.30 for feature ∆ (subscripts A and D are used to denote PMA and dose, resp.). Since PMA and dose are correlated (Pearson correlation coefficient rAD=0.45), we also define the coefficient of partial determination in order to account for the effect of both predictor variables on features S and ∆. Numerical values are equal to RA∣D2=0.002 and RD∣A2=0.49 for feature S and RA∣D2=0.05 and RD∣A2=0.20 for feature ∆. The statistical significance of the coefficients of (partial) determination is denoted in Figure 8. Finally, it is important to note that PMA and dose are not collinear using a linear model. This can be assessed by computing the variance inflatable factor (VIF) [35], which is equal to VIF=1.2572. A VIF close to 1 indicates the lack of collinearity.Figure 6
Signal interaction between MABP and EEG after administration of propofol att = 0 minutes. The signal interaction was computed using kTxi,xj. A reduction in interaction is observed among the different signal modalities after the administration of propofol, with a slow recovery to the reference level. The black line and gray shaded area present the median and interquartile range (IQR), respectively.Figure 7
Changes in the physiological network, assessed using a graph modelGT, for one neonate in the group of study (PMA 30 weeks, 1.0 mg·kg−1) at three different time instances: plots (a), (b), and (c) illustrate the edge weights for t = 10, 30, and 180 minutes after propofol administration, respectively. The graph model was constructed using kTxi,xj defined in (2). Under each graph, the average graph degree δGT is presented in a time frame starting right after propofol administration (t = 0) up to 6 hours after. The average graph degree measures the average connection strength of the graph edges. From (a) to (c), the edge weights increase, which translates in an increased δGT.
(a)
10 minutes after propofol
(b)
30 minutes after propofol
(c)
3 hours after propofolFigure 8
The relation between featuresS and ∆, computed from the MABP-EEG interaction curves presented in Figure 6, and PMA and propofol. The data points and the linear least squares fit are depicted in black and gray, respectively. The shaded area indicates the 95-percentage confidence bounds on the least squares fit. The coefficient of (partial) determination is indicated in each plot (subscripts A and D denote PMA and propofol dose, resp.). A single asterisk, double asterisks, and triple asterisks denote a p value smaller than 0.05, 0.01, and 0.001, respectively.
(a)
Relation betweenS and PMA
(b)
Relation betweenS and dose
(c)
Relation between ∆ and PMA
(d)
Relation between ∆ and dose


### 4.2. Overall Interactions

Figure9(a) presents a comparison of the vertex degree di (in red is the interaction of modality i with respect to the other modalities) with the average degree δGT (in black is the average interaction of all signal modalities) for all of the signal modalities after administration of propofol at t = 0 minutes. The curves are computed from graph models constructed using the similarity measure kTxi,xj (2). The results are presented for the whole group of study (N=22). Propofol-induced loss of interaction among the signals is associated with a drop in δGT. The drop in average graph degree can also be observed in Figure 7, which illustrates the graph model for one neonate in the group of study at different time instances. As shown in Figure 9(a), the δGT value is highly determined by dMABP during the first 30 minutes. Indeed, the MABP vertex degree is considerably lower compared to the degree of the other modalities in this time frame. From 30 minutes onwards, the increase of δGT to the reference level is highly influenced by dEEG, which is associated with the slowest recovery in signal dynamics.Figure 9
Comparing the vertex degree values (red) with the graph average degree (black) after administration of propofol att = 0 minutes. The graph models were constructed using kTxi,xj (a) and kFxi,xj (b). The results are presented for the whole group of study (N=22). From top to bottom, the vertex degree di is compared to the graph average degree δG for HR, MABP, SaO2, rScO2, and EEG, respectively. dMABP highly determines the signal interaction pattern during the first 30 minutes after propofol administration, while dEEG highly influences the signal interaction pattern from 30 minutes to 90 minutes after propofol administration. After 90 minutes, the neonates are recovered from propofol, as indicated by the steady reference levels observed after 90 minutes.
(a)
Graphs constructed usingkTxi,xj
(b)
Graphs constructed usingkFxi,xjFigure9(b) shows the vertex degree di (red) with the graph average degree δGF (black) after propofol administration at t = 0 for the graph models constructed using the second similarity measure, that is, kFxi,xj (3). As before, the results are presented for the whole group of study (N=22). A reduction in interaction can be observed after propofol administration, which is in agreement with the results of Figure 9(a). Again, MABP is observed to be the contributing factor in the propofol-induced loss of interaction during the first 30 minutes after propofol administration. Indeed, this vertex is associated with the lowest degree values during this time frame. From 30 minutes onwards, the increase of δGF is again influenced by EEG dynamics. This effect is however less pronounced compared to the observation of Figure 9(a). In general, the results from kTxi,xj and kFxi,xj are similar, which might indicate that time delayed and/or interaction of opposite signs are not present in our dataset or that the influence of those interactions is not relevant, probably due to the length of the analysis window (15 minutes) that we used in the analysis.


## 4.1. MABP-EEG Pairwise Interaction

The interaction curves of MABP with respect to EEG after administration of propofol att = 0 minutes are illustrated in Figure 6. These curves are computed using kTxi,xj, defined in (2). The EEG signal is represented by the running IBI duration, as outlined in Section 3. From top to bottom, the interaction pattern is shown for the entire group of study (N=22) and the individual age groups presented in Table 1. First, a pronounced loss in interaction is observed, followed by a gradual increase to a reference level, which is in general reached at t = 90 minutes. Note that this loss in interaction is present among all of the signal modalities of the multimodal dataset, as indicated by the graphs in Figure 7. Figure 8 presents the relation between the features used to quantify interaction strength (S and ∆) and PMA and propofol dose. In addition to the data points, the least squares linear fit is defined (straight lines), together with the 0.95 percentiles of the linear fit (shaded area). The goodness of the linear fit is assessed using the coefficient of determination Ri2, which is equal to RA2=0.09 and RD2=0.53 for feature S and RA2=0.17 and RD2=0.30 for feature ∆ (subscripts A and D are used to denote PMA and dose, resp.). Since PMA and dose are correlated (Pearson correlation coefficient rAD=0.45), we also define the coefficient of partial determination in order to account for the effect of both predictor variables on features S and ∆. Numerical values are equal to RA∣D2=0.002 and RD∣A2=0.49 for feature S and RA∣D2=0.05 and RD∣A2=0.20 for feature ∆. The statistical significance of the coefficients of (partial) determination is denoted in Figure 8. Finally, it is important to note that PMA and dose are not collinear using a linear model. This can be assessed by computing the variance inflatable factor (VIF) [35], which is equal to VIF=1.2572. A VIF close to 1 indicates the lack of collinearity.Figure 6
Signal interaction between MABP and EEG after administration of propofol att = 0 minutes. The signal interaction was computed using kTxi,xj. A reduction in interaction is observed among the different signal modalities after the administration of propofol, with a slow recovery to the reference level. The black line and gray shaded area present the median and interquartile range (IQR), respectively.Figure 7
Changes in the physiological network, assessed using a graph modelGT, for one neonate in the group of study (PMA 30 weeks, 1.0 mg·kg−1) at three different time instances: plots (a), (b), and (c) illustrate the edge weights for t = 10, 30, and 180 minutes after propofol administration, respectively. The graph model was constructed using kTxi,xj defined in (2). Under each graph, the average graph degree δGT is presented in a time frame starting right after propofol administration (t = 0) up to 6 hours after. The average graph degree measures the average connection strength of the graph edges. From (a) to (c), the edge weights increase, which translates in an increased δGT.
(a)
10 minutes after propofol
(b)
30 minutes after propofol
(c)
3 hours after propofolFigure 8
The relation between featuresS and ∆, computed from the MABP-EEG interaction curves presented in Figure 6, and PMA and propofol. The data points and the linear least squares fit are depicted in black and gray, respectively. The shaded area indicates the 95-percentage confidence bounds on the least squares fit. The coefficient of (partial) determination is indicated in each plot (subscripts A and D denote PMA and propofol dose, resp.). A single asterisk, double asterisks, and triple asterisks denote a p value smaller than 0.05, 0.01, and 0.001, respectively.
(a)
Relation betweenS and PMA
(b)
Relation betweenS and dose
(c)
Relation between ∆ and PMA
(d)
Relation between ∆ and dose


## 4.2. Overall Interactions

Figure9(a) presents a comparison of the vertex degree di (in red is the interaction of modality i with respect to the other modalities) with the average degree δGT (in black is the average interaction of all signal modalities) for all of the signal modalities after administration of propofol at t = 0 minutes. The curves are computed from graph models constructed using the similarity measure kTxi,xj (2). The results are presented for the whole group of study (N=22). Propofol-induced loss of interaction among the signals is associated with a drop in δGT. The drop in average graph degree can also be observed in Figure 7, which illustrates the graph model for one neonate in the group of study at different time instances. As shown in Figure 9(a), the δGT value is highly determined by dMABP during the first 30 minutes. Indeed, the MABP vertex degree is considerably lower compared to the degree of the other modalities in this time frame. From 30 minutes onwards, the increase of δGT to the reference level is highly influenced by dEEG, which is associated with the slowest recovery in signal dynamics.Figure 9
Comparing the vertex degree values (red) with the graph average degree (black) after administration of propofol att = 0 minutes. The graph models were constructed using kTxi,xj (a) and kFxi,xj (b). The results are presented for the whole group of study (N=22). From top to bottom, the vertex degree di is compared to the graph average degree δG for HR, MABP, SaO2, rScO2, and EEG, respectively. dMABP highly determines the signal interaction pattern during the first 30 minutes after propofol administration, while dEEG highly influences the signal interaction pattern from 30 minutes to 90 minutes after propofol administration. After 90 minutes, the neonates are recovered from propofol, as indicated by the steady reference levels observed after 90 minutes.
(a)
Graphs constructed usingkTxi,xj
(b)
Graphs constructed usingkFxi,xjFigure9(b) shows the vertex degree di (red) with the graph average degree δGF (black) after propofol administration at t = 0 for the graph models constructed using the second similarity measure, that is, kFxi,xj (3). As before, the results are presented for the whole group of study (N=22). A reduction in interaction can be observed after propofol administration, which is in agreement with the results of Figure 9(a). Again, MABP is observed to be the contributing factor in the propofol-induced loss of interaction during the first 30 minutes after propofol administration. Indeed, this vertex is associated with the lowest degree values during this time frame. From 30 minutes onwards, the increase of δGF is again influenced by EEG dynamics. This effect is however less pronounced compared to the observation of Figure 9(a). In general, the results from kTxi,xj and kFxi,xj are similar, which might indicate that time delayed and/or interaction of opposite signs are not present in our dataset or that the influence of those interactions is not relevant, probably due to the length of the analysis window (15 minutes) that we used in the analysis.


## 5. Discussion

In the present analysis, we study how different physiologic systems dynamically interact and collectively behave after a propofol bolus administration in preterm neonates. These physiologic systems are presented by the different signal modalities under study. Note that we focus on the interaction between the brain and the cardiovascular system. This study can therefore be situated in the interdisciplinary field of network physiology [5].Results indicate that propofol causes a change in the dynamical interactions between the different signals up to 90 minutes after propofol administration. The strength of this effect was observed to be mainly determined by propofol dose. In addition, the recovery phase was observed to be mainly determined by EEG dynamics, due to a much slower recovery to the reference level compared to the other signal modalities.


### 5.1. MABP-EEG Pairwise Interaction

Sedation of neonates using propofol induces a reduction in the interaction between MABP and EEG (Figure6), with only a slow, gradual increase back to the reference level. The most pronounced decrease in interaction pattern is associated with the oldest neonates in the group of study (moderate to late preterm): a strong loss of interaction is observed during the first 60 minutes after propofol administration, followed by a brisk increase back to baseline (Figure 6(d)). This pattern clearly differs from that of the younger neonates (extremely to very preterm), which are characterized by a less-pronounced reduction in interaction and a more gradual increase back to reference levels (Figures 6(b) and 6(c)).Two possible indicators for the observed difference in signal interaction patterns are proposed. Both indicators are based on signal amplitude changes, since the signal interaction measurekTxi,xj highly depends on signal amplitudes. Firstly, the discontinuity pattern of neonatal EEG changes with age. Especially, the oldest neonates (moderate to late preterm) are characterized by a much more continuous EEG pattern (tracé continue) compared to the younger neonates (extremely to very preterm; tracé discontinue) [30]. A more continuous EEG can result in a more pronounced increase in IBI duration after propofol, potentially explaining the more pronounced loss in signal interaction observed among the oldest neonates in the group of study. Secondly, Simons et al. observed a higher incidence of hypotension with increasing dose of propofol [10]. In this study, higher doses were administered to older neonates, as demonstrated by Table 1. Evidently, a more pronounced impact on MABP can be responsible for a stronger loss in signal interaction.Since PMA and propofol dose (predictor variables) are correlated (rAD=0.45), the influence of each factor on the resulting signal interaction pattern is assessed using features S and ∆ (independent variables). Figure 8 presents the relation between these features and PMA and propofol dose. From Figure 8, it is clear that the influence of PMA on the independent variables is minimal, especially when taking into account the influence of the dose. Indeed, the coefficients of partial determination are very small for PMA. (RA∣D2=0.002 and RA∣D2=0.05 for S and ∆, resp.). This observation is confirmed by the fact that the coefficient of partial determination is only slightly smaller compared to the coefficient of determination for propofol dose, especially for feature S. Therefore, it is clear that the interaction between MABP and EEG is mainly influenced by propofol dose. The difference in interaction pattern observed in Figure 6 is thus mainly caused by the difference in propofol dose administered to the neonates in the different age groups, and not by the difference in PMA.


### 5.2. Overall Interactions

The phase of sedation using propofol is characterized by a markedly different network structure compared to the reference phase, indicating a clear association between network topology and physiologic function. This is illustrated in Figure7: after 10 minutes, the graph is weakly connected indicating a highly reduced overall signal interaction as opposed to the strongly connected graph observed at 3 hours after propofol administration.MABP is observed to be the main contributor to the reduction in signal interaction during the first 30 minutes after propofol administration, as indicated in Figure9. During this time frame, MABP strongly influences the strength of the overall interaction pattern, since the vertex degree is lower compared to the average graph degree. This effect can partly be explained as an amplitude effect. Indeed, propofol administration is associated with a pronounced decrease in MABP, which can last up to one hour after propofol administration, as described by many authors [6–8, 10]. The physiologic response of the other signal modalities is less affected by propofol compared to MABP. This pronounced change in signal amplitude could explain why MABP highly influences the overall interactions, especially during the first 30 minutes after propofol administration. It is important to note, however, that the explained loss in signal interaction can not be entirely explained by only taking into account the signal amplitude and change thereof in time. Indeed, the propofol-induced loss in signal interaction is also observed in Figure 9(b), which presents the results using similarity measure kFxi,xj. This measure assesses the interaction of the signals in the frequency domain.For 30 minutes up to 90 minutes after propofol administration, the degree of the EEG signal is considerably lower than the degree values of the other modalities. As before, this finding can be observed in Figure9. The EEG signal is the only signal associated with degree values below the average degree, indicating the slow recovery of EEG dynamics with respect to the other modalities. Thus, MABP dynamics recover faster (generally recovered 30 minutes after propofol administration) compared to EEG dynamics (recovery takes up to 90 minutes after propofol administration). From a signal processing point of view, this might indicate the safety of propofol, since MABP can adapt to the needs of brain metabolism, once the EEG signal is recovered. It is important to note, however, that the neonates included in the present analysis were all sedated using propofol as part of an INSURE procedure. Surfactant causes a significant decrease in EEG activity, which can last up to 24 hours after surfactant administration, as described by van den Berg et al. [36]. Therefore, surfactant could also influence the decreased EEG interactions observed in Figure 9. The extent of this effect is however not clear at this point, since no control group without surfactant was available to compare with.From 90 minutes after propofol administration onwards, the vertex degree and average degree curves presented in Figure9 are characterized by stable reference levels. This indicates that the signal interaction pattern is restored after propofol administration.


## 5.1. MABP-EEG Pairwise Interaction

Sedation of neonates using propofol induces a reduction in the interaction between MABP and EEG (Figure6), with only a slow, gradual increase back to the reference level. The most pronounced decrease in interaction pattern is associated with the oldest neonates in the group of study (moderate to late preterm): a strong loss of interaction is observed during the first 60 minutes after propofol administration, followed by a brisk increase back to baseline (Figure 6(d)). This pattern clearly differs from that of the younger neonates (extremely to very preterm), which are characterized by a less-pronounced reduction in interaction and a more gradual increase back to reference levels (Figures 6(b) and 6(c)).Two possible indicators for the observed difference in signal interaction patterns are proposed. Both indicators are based on signal amplitude changes, since the signal interaction measurekTxi,xj highly depends on signal amplitudes. Firstly, the discontinuity pattern of neonatal EEG changes with age. Especially, the oldest neonates (moderate to late preterm) are characterized by a much more continuous EEG pattern (tracé continue) compared to the younger neonates (extremely to very preterm; tracé discontinue) [30]. A more continuous EEG can result in a more pronounced increase in IBI duration after propofol, potentially explaining the more pronounced loss in signal interaction observed among the oldest neonates in the group of study. Secondly, Simons et al. observed a higher incidence of hypotension with increasing dose of propofol [10]. In this study, higher doses were administered to older neonates, as demonstrated by Table 1. Evidently, a more pronounced impact on MABP can be responsible for a stronger loss in signal interaction.Since PMA and propofol dose (predictor variables) are correlated (rAD=0.45), the influence of each factor on the resulting signal interaction pattern is assessed using features S and ∆ (independent variables). Figure 8 presents the relation between these features and PMA and propofol dose. From Figure 8, it is clear that the influence of PMA on the independent variables is minimal, especially when taking into account the influence of the dose. Indeed, the coefficients of partial determination are very small for PMA. (RA∣D2=0.002 and RA∣D2=0.05 for S and ∆, resp.). This observation is confirmed by the fact that the coefficient of partial determination is only slightly smaller compared to the coefficient of determination for propofol dose, especially for feature S. Therefore, it is clear that the interaction between MABP and EEG is mainly influenced by propofol dose. The difference in interaction pattern observed in Figure 6 is thus mainly caused by the difference in propofol dose administered to the neonates in the different age groups, and not by the difference in PMA.


## 5.2. Overall Interactions

The phase of sedation using propofol is characterized by a markedly different network structure compared to the reference phase, indicating a clear association between network topology and physiologic function. This is illustrated in Figure7: after 10 minutes, the graph is weakly connected indicating a highly reduced overall signal interaction as opposed to the strongly connected graph observed at 3 hours after propofol administration.MABP is observed to be the main contributor to the reduction in signal interaction during the first 30 minutes after propofol administration, as indicated in Figure9. During this time frame, MABP strongly influences the strength of the overall interaction pattern, since the vertex degree is lower compared to the average graph degree. This effect can partly be explained as an amplitude effect. Indeed, propofol administration is associated with a pronounced decrease in MABP, which can last up to one hour after propofol administration, as described by many authors [6–8, 10]. The physiologic response of the other signal modalities is less affected by propofol compared to MABP. This pronounced change in signal amplitude could explain why MABP highly influences the overall interactions, especially during the first 30 minutes after propofol administration. It is important to note, however, that the explained loss in signal interaction can not be entirely explained by only taking into account the signal amplitude and change thereof in time. Indeed, the propofol-induced loss in signal interaction is also observed in Figure 9(b), which presents the results using similarity measure kFxi,xj. This measure assesses the interaction of the signals in the frequency domain.For 30 minutes up to 90 minutes after propofol administration, the degree of the EEG signal is considerably lower than the degree values of the other modalities. As before, this finding can be observed in Figure9. The EEG signal is the only signal associated with degree values below the average degree, indicating the slow recovery of EEG dynamics with respect to the other modalities. Thus, MABP dynamics recover faster (generally recovered 30 minutes after propofol administration) compared to EEG dynamics (recovery takes up to 90 minutes after propofol administration). From a signal processing point of view, this might indicate the safety of propofol, since MABP can adapt to the needs of brain metabolism, once the EEG signal is recovered. It is important to note, however, that the neonates included in the present analysis were all sedated using propofol as part of an INSURE procedure. Surfactant causes a significant decrease in EEG activity, which can last up to 24 hours after surfactant administration, as described by van den Berg et al. [36]. Therefore, surfactant could also influence the decreased EEG interactions observed in Figure 9. The extent of this effect is however not clear at this point, since no control group without surfactant was available to compare with.From 90 minutes after propofol administration onwards, the vertex degree and average degree curves presented in Figure9 are characterized by stable reference levels. This indicates that the signal interaction pattern is restored after propofol administration.


## 6. Conclusions

In this study, we have shown that graph theory can be used to assess changes in signal interaction and that the resulting graph models can be used to study the difference between distinct physiologic states.Moreover, for our propofol case study, we derived that the overall signal interaction pattern after propofol administration is highly influenced by both MABP and EEG. The MABP signal is the main contributor to the loss in signal interactions during the first 30 minutes after propofol, due to the strong decoupling of MABP dynamics with respect to the other signal modalities, while the EEG signal highly influences the interaction pattern thereafter. This finding indicates that MABP dynamics recover first, followed by a much slower recovery of the EEG signal, meaning that MABP dynamics are recovered while EEG metabolism is still down. Thus, when EEG dynamics recover, MABP can adapt to supply new needs of the brain in order to sustain its function.Propofol affects signal dynamics with an overall recovery time of around 90 minutes, as assessed by the graph average degree. After 90 minutes, these curves are characterized by steady reference levels, indicating that, at least from a biosignal processing point of view, the overall signal dynamics are recovered from propofol and that the physiological system is associated with a high degree of signal interaction.The signal interaction pattern observed after propofol administration is influenced only by propofol dose, and thus not by PMA. This relation was observed for the pairwise interaction curves and the system interaction measure (average graph degree) derived from the graph model of the neonate.



---

*Source: 6504039-2018-08-09.xml*