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

Department of Electrical Engineering (ESAT), STADIUS, KU Leuven, Leuven, Belgium imec, Leuven, Belgium Department of Development and Regeneration, KU Leuven, Leuven, Belgium Department of Neonatology, UZ Leuven, Leuven, Belgium Department of Pediatric Surgery and Intensive Care, Erasmus MC-Sophia Children’s Hospital, Rotterdam, Netherlands Department of Neonatology, Erasmus MC-Sophia Children’s Hospital, Rotterdam, Netherlands


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 (rScO 2 ), 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 rScO 2 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 (SaO 2 ), rScO 2 , 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][7][8][9][10][11], children [12], and adults [13][14][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, SaO 2 , and rScO 2 , 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. Section 2 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.

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)(27)(28)(29)(30)(31)(32)(33)(34)(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   2 Complexity about the clinical characteristics of the subjects can be found in the original paper by Smits et al. [6]. Practices 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][10][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, SaO 2 , rScO 2 , 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 threecompartment 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.
The systemic variables (HR (beats/min), MABP (mmHg), and SaO 2 (%)) 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 rScO 2 (%) noninvasively with INVOS 5100 using a cerebral neonatal OxyAlert NIRS sensor (Covidien, Mansfield, Massachusetts). As for the systemic variables, the sampling frequency for rScO 2 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.

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 3 Complexity 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.

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 Figure 3. 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, A complete graph with n vertices has m = n n − 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,  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 k T x i , x j uses the raw signals in the RBF kernel and is thus defined as where x i and x j 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, x i and x j are segments with a length of 15 minutes, as mentioned before. The similarity k T x i , x j is bounded by 0 (absence of common interactions) and 1 (exact   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. 5 Complexity 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 k T x i , x j is denoted as G T .
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 measure k F x i , x j is defined as where S x i and S x j represent the PSD of input signals x i and x j (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 k F x i , x j 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 k T x i , x j 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 k F x i , x j is denoted as G F .

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 k T x i , x j and k F x i , x j both depend on this parameter σ. The optimization procedure is the same for both similarity measures. Therefore, it is outlined in terms of k x i , x j , 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 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 where p k 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 : with 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). Figure 4 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)(27)(28)(29) and 94 (86-113) for k T x i , x j and k F x i , x j , respectively.

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 matrix A of a graph, the degree d i of a vertex, and the average degree δ G of a graph G.

Adjacency Matrix.
A weighted graph G consists of a nonempty finite set V of elements called vertices v i (or nodes) and a finite set E of distinct unordered pairs of distinct elements of V called edges w ij (or links) [33]. Note that the edges of the graph are represented by their weights w ij . 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 v i of V and their corresponding 6 Complexity weights. More precisely, the adjacency matrix A is constructed as if there is an edge between υ i and υ j , 0, otherwise 8

Vertex
Degree. The degree d j associated with a vertex v j of an undirected weighted graph G, with adjacency matrix A, is defined as the sum of all edges incident to v j : where n is the number of vertices. Therefore, the degree d j characterizes the connection strength of the vertex v j 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 w ij ∈ 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, If d j = 0, vertex v j is called an isolated vertex, since it is not connected to any other vertex of the graph. A vertex degree and 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 d j 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.   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.

Complexity
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.
In the present analysis, the relation between the feature values S and Δ (dependent variables) and PMA and propofol dose (predictor variables) is studied using linear regression models. The coefficient of determination R i 2 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 10 5 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 Math-Works, Natick, Massachusetts). Graph theory analysis is performed using the MATLAB toolbox for network analysis, provided by MIT Strategic Engineering [34].

MABP-EEG Pairwise
Interaction. The interaction curves of MABP with respect to EEG after administration of propofol at t = 0 minutes are illustrated in Figure 6. These curves are computed using k T x i , x j , 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 R 2 i , which is equal to R 2 A = 0 09 and R 2 D = 0 53 for feature S and R 2 A = 0 17 and R 2 D = 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 r AD = 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 R 2 A|D = 0 002 and R 2 D|A = 0 49 for feature S and R 2 A|D = 0 05 and R 2 D|A = 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 V IF = 1 2572. A VIF close to 1 indicates the lack of collinearity. Figure 9(a) presents a comparison of the vertex degree d i (in red is the interaction of modality

Complexity
i with respect to the other modalities) with the average degree δ G T (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 k T x i , x j (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 δ G T . 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 δ G T value is highly determined by d MABP 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 δ G T to the reference level is highly influenced by d EEG , which is associated with the slowest recovery in signal dynamics. Figure 9(b) shows the vertex degree d i (red) with the graph average degree δ G F (black) after propofol administration at t = 0 for the graph models constructed using the second similarity measure, that is, k F x i , x j (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.  Figure 6: Signal interaction between MABP and EEG after administration of propofol at t = 0 minutes. The signal interaction was computed using k T x i , x j . 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.

Complexity
Indeed, this vertex is associated with the lowest degree values during this time frame. From 30 minutes onwards, the increase of δ G F 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 k T x i , x j and k F x i , x j 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.

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.

MABP-EEG Pairwise
Interaction. Sedation of neonates using propofol induces a reduction in the interaction between MABP and EEG ( Figure 6), 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 measure k T x i , x j 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 (r AD = 0 45), the influence of each factor on the resulting signal interaction pattern is assessed using features  10 Complexity 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. (R 2 A|D = 0 002 and R 2 A|D = 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.

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 Figure 7: after 10 minutes, the graph is   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. 11 Complexity 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 Figure 9. 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][7][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 k F x i , x j . 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 Figure 9. The EEG signal is the only signal associated with degree values below the average : Comparing the vertex degree values (red) with the graph average degree (black) after administration of propofol at t = 0 minutes. The graph models were constructed using k T x i , x j (a) and k F x i , x j (b). The results are presented for the whole group of study (N = 22). From top to bottom, the vertex degree d i is compared to the graph average degree δ G for HR, MABP, SaO 2 , rScO 2 , and EEG, respectively. d MABP highly determines the signal interaction pattern during the first 30 minutes after propofol administration, while d EEG 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. 12 Complexity 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 Figure 9 are characterized by stable reference levels. This indicates that the signal interaction pattern is restored after propofol administration.

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.

Data Availability
The data used to support the findings of this study are restricted by the Ethische Commissie onderzoek UZ/KU Leuven in order to protect patient privacy. Data are available from Dries Hendrikx (dries.hendrikx@esat.kuleuven.be) for researchers who meet the criteria for access to this confidential data.

Disclosure
This paper reflects only the authors' views and the union is not liable for any use that may be made of the contained information.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper, since the received funding, as stated in the Acknowledgments, does not lead to any conflicts of interest.