Deoxyhemoglobin

We quantified the anatomical connectivity using graph theoretical measures [23] where the in-degree and out-degree are the number of incoming and outgoing connections to/from a node. The degree is the sum of in- and out-degree. The clustering coefficient is the number of all existing connections between a node's neighbors divided by all such possible connections. The betweeness centrality is the fraction of the shortest path between any two pairs of nodes passing through a particular node.
The network model with the coupling term of strength c is implemented as: where u_i, ν_i are the state variables of the ith neural population and f_ij is the connectivity matrix. White Gaussian noise n_u(t), n_ν(t) is introduced additively. The functions g and h are based on FitzHugh-Nagumo systems [26] (link),[27] with and h(u_i,ν_i) = −(1/τ)[u_i−α+bν_i], and α = 1.05, β = 0.2, γ = 1.0, τ = 1.25. For the stability analysis (no noise) we employed Matlab DDE23 to solve the coupled delay differential equations. The coupled delay differential equations with additive noise were solved in Matlab by a simplified and faster algorithm. More specifically, we employed a standard fourth order Runge-Kutta method for integrating the intrinsic Fitz-Hugh Nagumo dynamics while the coupling and the stochastic terms were integrated using Euler method. The step size for the simulation was 0.001 and we confirmed that no better convergence of solution was achieved using smaller step sizes to ensure numerical convergence.
The time delays are computed from the Euclidean distance matrix d_ij of the locations of the brain areas i and j. To do so, the three-dimensional regional map locations were converted to approximate Talairach stereotaxic atlas locations by first identifying the mapping of regional map locations as designated on the human brain to the anatomical locations in Talairach space using the Anatomical Automatic Labeling (AAL) image provided by Tzourio-Mazoyer et al. [36] (link). Once the approximate location was identified in the AAL brain, the coordinate for the centre of the AAL region was used for the location of the corresponding regional map location. Each region was represented as a surface composed of a sufficient number of triangles. To obtain the triangulation, a T1-weighted MR image from a single human subject was segmented in grey and white matter compartments and the cortical surface represented as a triangular net using the CURRY software package (Compumedics Neuroscan, Ltd). The T1 image was co-registered to a standard MRI atlas (MNI305, [37] ) using a 12-parameter affine transform with sinc interpolation as implemented in SPM99 (see http://www.fil.ion.ucl.ac.uk/spm/ and [38] ). The transform matrix from the co-registration was then applied to the triangulated cortical surface to the MRI atlas.
The stability diagram for the network in Equation (1) is obtained by linear stability analysis leading to the characteristic equation where ū is the fixed point solution. The eigenvalue λ has 2N non-trivial roots with The equilibrium state is stable if all eigenvalues λ have negative real parts, Re(λ)<0, which were found numerically. The stability diagrams in Figure 3 were constructed using this procedure. We also cross-validated the presence of negative real parts of the eigenvalues by direct numerical simulations of Equation 1.
We obtain activity at different areas by simulating Equation 1 for parameter values indicated in the stability diagram. The parameters are chosen to lie on or just below the critical boundary of stable and unstable regions. Network data are simulated for numerical values of parameters in the stable region. Once the network dynamics settles into its equilibrium state (see Figure S1), the coupling parameter c is increased just beyond the critical boundary. We use the smallest increase of c possible given the discretization of the parameter space. As a consequence, now in the unstable regime, the network dynamics increases towards high-amplitude oscillations. A typical time series plot is shown in Figure S1. Using a sliding temporal window of 500 ms width, we perform a Principal Component Analysis (PCA) during the transient as the oscillations increase. The local Center Manifold Theorem guarantees that the network modes with the largest positive real part of the eigenvalue grow fastest and hence dominate the transient initially. Hence, the eigenvectors of PCA span a linear vector space, in which the dominant network modes will be represented. In other words, the networks implicated in rest state activity will be a linear superposition of the PCA eigenvectors. Then the spatiotemporal data can be decomposed as: where the kth PCA eigenvector ψ_k spans a spatial network. During all transients observed in our simulations, the first two PCA eigenvectors contribute together at least 99.995 percent (see Figure 4B). Hence it does suffice to represent the entire transient dynamics by the first two PCA eigenvectors. Since, in a given PCA eigenvector ψ_k, each node is multiplied by the same time-dependent coefficient ξ_k(t), the magnitude of the ith vector element will scale the resulting contribution of the ith node to the network dynamics. The most dominant nodes of these two networks ψ_k(i) are then identified through an ordering process: we compute ψ_k(i)² for all nodes i and both network modes k = 1, 2 and order these according to power (see for instance Figure 4C). There is no hard criterion to identify a threshold for the inclusion of nodes in a network. For reasons of clarity, we choose to show the first three dominating nodes for each eigenvector in Figure 4A, which corresponds to at least 90% of the power per eigenvector in all cases.
To relate the simulated neural activity to recent fMRI studies, we have generated BOLD signal for each regions by using a hemodynamic model. This model combines the Balloon/Windkessel model comprised in venous volume and deoxyhemoglobin content with a linear dynamical model of how synaptic activity causes changes in regional cerebral blood flow [31] (link). For each region, neural activity causes an increase in a vasodilatory signal inducing blood flow, which changes blood volume and deoxyhemoglobin content. The BOLD signal is given by a volume-weighted sum of extra- and intra vascular signals as the function of volume and deoxyhemoglobin content. The local neural activity, which is taken to be the absolute value of the time derivative of the output occurring by our network model in each brain region, is used as the main model input to estimate a BOLD signal. For the analyses, the global mean signal (average over all regions) has been regressed out from the single BOLD time series. All parameters regarding blood flow, deoxyhemoglobin content, and vessel volume in the model equation are taken from [31] (link).

To conduct the signal separation using Eqs. (3) and (4), the coefficients of both hemodynamic modalities, and , have to be known. In the following section, we propose a procedure to determine them.
The functional component originates mainly from the regional cerebral hemodynamics evoked by neural activation. We surveyed fNIRS studies from the past and utilized studies satisfying two conditions: (1) those that used an appropriate experimental design to exclude physiological signals other than cerebral hemodynamics and (2) those that provided a graph of both oxy- and deoxyhemoglobin changes from the beginning to the end of functional activation. We calculated the value of from the graph. A vertically flipped copy of the graph was superposed on the original graph using a drawing software (Illustrator, Adobe Systems), and the original's vertical magnification was manually modulated such that the oxyhemoglobin changes in the original graph visually coincided with the deoxyhemoglobin changes in the copy. The magnification percentage gave . In most cases, the shape of the flipped deoxyhemoglobin graphs was quite similar to that of the original oxyhemoglobin graphs. This indicated that the negative linear correlation between oxy- and deoxyhemoglobin generally held true for the hemodynamics associated with neural activation. The result is shown in Table 1. All values given by the studies fell in a range of . In Table 1, the mean SD of was −0.56 0.12. From the statistical analysis using ANOVA, the values in Table 1 showed a statistical difference between groups with visual and other kinds of stimulation ( for visual and for other kinds of stimulation; ). There was no statistical difference in between different subjects and durations. From these analyses, we adopted as the universal value.
The coefficient of the systemic component has the relationship with . It is known that oxygen saturation levels differ among vessels. For example, under normal conditions, the saturation level in arteries is greater than 95%, while that in veins is approximately 70% [46] (link), and it decreases to approximately 50% after intense physical activity. These saturation levels correspond to of less than 0.053, 0.43, and 1.0, respectively.
The coefficient may vary according to the type of task. If a task containing a psychophysiological load is executed, it may change the arterial blood pressure, respiration rate, and vasomotor action. These changes cause blood volume changes mainly in the arteries and arterioles. Since blood in the arteries and arterioles has a high oxygen saturation level, should be small. In contrast, hyperemia induced by a posture change gives a larger value because a passive volume capacity change can occur not only in arteries but also in veins having lower oxygen saturation levels. Furthermore, different values are expected if the intensity of physical activity is varied. Thus, unlike , we cannot expect a universal value of under various task conditions.
In our model, we assume that capillaries mainly generate the functional component and arteries and veins generate the systemic component. Since these components originate from different vessels and different hemodynamic modalities, we assumed a high statistical independence between them. Thus, we determined by minimizing the mutual information between these components. The mutual information is given as follows. where and represent the probability density functions of and , respectively, and represents a joint probability density function of and . These probabilities are estimated by the normalized histograms of and , and the normalized joint histograms of and . Histograms are calculated from Eqs. (3) and (4) when and are given. Here we fix based on the discussion on Table 1. Thus, if we set , the mutual information can be calculated by Eq. (5). By enumerating in we determine , which minimizes the mutual information.

Figures 8(a)–8(d) show the evaluation results of the smartwatch-based prototype according to different types of oral challenges for the same male subject who has no underlying health conditions. Here, the vertical dashed line in each plot indicates the time when the oral challenge was given to the subject, and the 10% vertical error bar range was applied to the reference CGM values considering the typical accuracy of the sensor.³⁴ (link)
In Fig. 8(a), the metabolic index MI agreed well with the reference glucose concentration transition induced by the sugary drink ingestion. On the other hand, in Figs. 8(c) and 8(d), both MI and reference values were maintained nearly flat during each oral challenge test. Obviously, the sugar-free drink does not affect the subject’s BGL that much. However, the finding that MI values did not change as much during the sugar-free oral challenge tests suggests that MI is not simply a reflection of the body’s hydration status. At the same time, judging from the results that both carbonated and non-carbonated sugar-free drinks showed similar flat MI trends, the result in Fig. 8(a) is unlikely to be due to the carbonation status of the oral challenge. Regarding Fig. 8(b), an offset of about 15 min that can be seen between MI and CGM is likely due to the lag between blood glucose and interstitial fluid (ISF) glucose. In general, 15 min of CGM sensor delay is a bit large compared to the typical delay reported by the manufacturer. However, a previous study³⁵ (link) shows that CGM delay can reach 15 min. Then, by compensating for 15 min of offset in Fig. 8(b), the plot can be modified as in Fig. 9(a). In this case, the MI also follows the up and down trend of the CGM well, in the same way as seen in Fig. 8(a). The proposed MI method response time can be estimated as 5 to 15 min from Figs. 8(a)–8(d) considering that the CGM device time delay is 5 to 15 min and the MI exhibited a similar trend.
By compiling the plot data from Figs. 8(a), 8(c), 8(d), and 9(a), a scatter plot Fig. 9(b) can be obtained. Here, each numerical data of MI corresponding to each CGM data point has been calculated by interpolation. In Fig. 9(b), the correlation coefficient $r$ calculated by the linear least squares (LLS) fitting reached 0.78, which indicates that the proposed metabolic index MI has a strong correlation with BGL.
In addition, the arterial oxygen saturation ${\mathrm{SaO}}_2(t)$ remained between $90\pm 2\%$ during a series of experiments, which means that ${\mathrm{SaO}}_2(t)·[1-{\mathrm{SaO}}_2(t)]$ remained between 0.07 and 0.11 in Eq. (32), and the maximum to minimum ratio of ${\mathrm{SaO}}_2(t)·[1-{\mathrm{SaO}}_2(t)]$ was at most $\sim 1.4$ . Here, the decrease in ${\mathrm{SaO}}_2(t)$ from its normal range is attributed to the measurement section where ${\mathrm{SaO}}_2(t)$ approaches ${\mathrm{StO}}_2(t)$ .³⁶ (link)
Therefore, the increase in the metabolic index MI was largely due to the change in $\mathrm{\Delta }\theta (t)$ : the phase delay between oxy- and deoxyhemoglobin. For reference, the $\mathrm{\Delta }\theta (t)$ transition in the experiment of Fig. 8(a) is shown in Fig. 9(c). In this study, the results consistently showed that $\mathrm{\Delta }\theta (t)$ was positive during each oral challenge test, indicating that the oxyhemoglobin NIRS signal $N_{{\mathrm{HbO}}_2}(t)$ consistently precedes the deoxyhemoglobin NIRS signal $N_{\mathrm{Hb}}(t)$ . The MI trend is therefore almost identical to the non-absolute $\mathrm{\Delta }\theta (t)$ trend. For this reason, only the MI trends are plotted in this paper to evaluate the proposed method, rather than separately plotting $\mathrm{\Delta }\theta (t)$ and ${\mathrm{SaO}}_2(t)$ .
Thus, the validity of the assumptions in Eq. (24) is verified.