# Numerical Model for Cerebrovascular Hemodynamics with Indocyanine Green Fluorescence Videoangiography

## Article information

## Abstract

### Objective

The use of indocyanine green videoangiography (ICG-VA) to assess blood flow in the brain during cerebrovascular surgery has been increasing. Clinical studies on ICG-VA have predominantly focused on qualitative analysis. However, quantitative analysis numerical modelling for time profiling enables a more accurate evaluation of blood flow kinetics. In this study, we established a multiple exponential modified Gaussian (multi-EMG) model for quantitative ICG-VA to understand accurately the status of cerebral hemodynamics.

### Methods

We obtained clinical data of cerebral blood flow acquired the quantitative analysis ICG-VA during cerebrovascular surgery. Varied asymmetric peak functions were compared to find the most matching function form with clinical data by using a nonlinear regression algorithm. To verify the result of the nonlinear regression, the mode function was applied to various types of data.

### Results

The proposed multi-EMG model is well fitted to the clinical data. Because the primary parameters—growth and decay rates, and peak center and heights—of the model are characteristics of model function, they provide accurate reference values for assessing cerebral hemodynamics in various conditions. In addition, the primary parameters can be estimated on the curves with partially missed data. The accuracy of the model estimation was verified by a repeated curve fitting method using manipulation of missing data.

### Conclusion

The multi-EMG model can possibly serve as a universal model for cerebral hemodynamics in a comparison with other asymmetric peak functions. According to the results, the model can be helpful for clinical research assessment of cerebrovascular hemodynamics in a clinical setting.

**Keywords:**Cerebrovascular circulation; Vascular surgical procedures; Indocyanine green; Fluorescence angiography; Hemodynamic monitoring; Computer-assisted numerical analysis

## INTRODUCTION

Indocyanine green videoangiography (ICG-VA) is a renowned tool for providing relevant information on vascular blood flow of patients during cerebrovascular surgery [5,10,13,18,29,31]. Several medical research groups have employed ICG-VA in neurosurgery studies. Most of them have used ICG-VA to assess vascular blood flow with qualitative analysis based on angiography findings. Only a few quantitative analyses on ICG-VA involved measurement of the fluorescence intensity in the process of time [1,4,19,27]. Tempolabile quantitative analysis enables more accurate assessments of hemodynamics and can lead to a model for blood flow in body [17,34]. Some studies have strived to provide several model functions for blood flow in the brain [3,9,22,24,35]. However, determining blood flow in a clinical setting is difficult because the signal for the vascular blood flow time profile is highly variable depending on the conditions of patients, as described in other imaging modalities [36]. Thus, it is important to establish a vascular blood flow model that can be applied to various time profile signals for data to elucidate the kinetics of blood circulation in the brain and obtain principle parameters that enable evaluation of biological phenomena. These parameters include the number, time (or location of the center), height, and growth and decay rates of the peaks [17]. These parameters provide objective information to diagnose the vascular blood flow of patients and identify the distinctions among the hemodynamic mechanisms [15,20,21,28].

In this study, we demonstrated that cerebrovascular hemodynamics can be presented as a mathematical model that is an exponential modified Gaussian (EMG) function. The numerical model was shown to be well fitted to ICG-VA tempolabile quantitative data. The model was constructed using a fitting algorithm with nonlinear regression. We predicted the primary parameters of curves (center of peaks, height of peaks, growth rate, and decay rate from the data set) of data comprised of points missed partly by the measurement limitation.

## MATERIALS AND METHODS

The study was approved by IRB of Seoul National University Boramae Hospital (30-2020-097). We obtained the curve data from ICG-VA during cerebrovascular surgery in conformity with our previous research [33]. We employed commercially available products for the ICG dye (Daiichi Sankyo, Tokyo, Japan) and surgical microscope (OPMI Pentero; Carl Zeiss Co., Oberkochen, Germany). The angiogram video was converted to static images using the video player and capture freeware (GomPlayer v. 2.1.36; Gretech Corp., Seoul, Korea). Tempolabile quantitative analysis data were extracted from the static image stacks by the multi-measure plug-in of ImageJ (v. 1.42q; http://rsb.info.nih.gov/ij/; National Institutes of Health, Bethesda, MD, USA).

The vascular blood flow circulatory system included the heart, lung, brain, and extra body [2]. ICG was injected as an intravenous bolus into the central venous line and was observed using an ICG-VA technique at the cerebral arteries. The measured signals were composed of absorption and elimination phases, and the time and width of the peaks in the time domain projected the kinetic information on the blood flow. Intravenous bolus injection formed an ideal pulse shape of the cerebrovascular hemodynamic curve between the absorption and elimination phases (Fig. 1A). However, the measured curve had the shape to an asymmetry peak function, as shown in Fig. 1B, because of the injection distribution and dispersion effect. The quantity of the injected ICG could have a Gaussian distribution in practice on account of the random noise variance. Thus, the ICG concentration in the blood could increase under the cumulative distribution function (CDF) in the absorption phase.

After the bolus, the concentration was exponentially decayed by the dispersion effect [14]. Therefore, we assumed that the measured cerebrovascular hemodynamic curve (CH_{measure}(t)) was represented as follow.

We chose the EMG model, which is one of the asymmetry peaks, as a candidate to fit on the curve. It has a similar mathematical express as Equation (1) [11]. The error function (erf) in the EMG could be transposed into CDF, Ø(t) with the following relationship.

EMG was compared and evaluated to identify a suitable model function for the measurement data. The data were calculated by the Levenberg-Marquardt algorithm (nonlinear least-squares algorithm) using scientific graphing and data analysis software, Origin9 Pro (v. 9.0.0.; OriginLab Co., Northampton, MA, USA).

### EMG model for the cerebrovascular hemodynamic curve

EMG has typically been used for both chromatographic distribution and peak analyses in a wide range of fields [8,12,23]. EMG is expressed by the convolution of one CDF Ø(t) and one exponential function (Equation [3]) :

where A is the peak amplitude, t_{0} is the exponential relaxation rate for decay, *w* is the variable related to the peak width and growth rate, and t_{C} is the variable related to the center of peaks.

The CDF with the sigmoid (S-shaped) form increases from 0 to 1. The parameter w is the growth steepness; thus, w means the growth rate degree. EMG is the convolution function. Meanwhile, CDF does not influence the elimination phase because the CDF value is only one in the elimination phase. In other words, CDF is dominant in the absorption phase, which is the growth part of the peak of the EMG curve. The exponential function is only related to the elimination phase, which is the decay part of the curve.

The decay rate degree enables defining the decay constant of the exponential function,

## RESULTS

To evaluate the validity of the EMG model, the R-squared value was compared with other peak functions and a bi-summation of the peak functions using the fitting algorithm. The applied data had a general form of a single peak, which was one of the quantitative ICG-VA clinical data. The number of frames in the X-axis of all measured data was transformed into the time domain to verify the data time profile.

We chose five peak functions : Gaussian function (G(t)), which is not actually an asymmetric peak function; nevertheless, it is the simplest function among the peak functions; log-normal function (L(t)) Gumbel distribution function (GB(t)); Maxwell-Boltzmann function (MB(t)); and EMG function (EMG(t)). Each function can be expressed as :

where y_{0}s are the offsets, A the amplitude of the peaks, w the width of the peaks, and t_{C} denotes the center of peaks.

### Comparison with other fitting functions

Fig. 2A shows the cerebrovascular hemodynamic data in the time domain. It is comprised of clinical data with a common form of a single peak in the quantitative ICG-VA with the time profile. Specifically, it consists of a peak and a long-range tail. We fitted the clinical data from the peak area (up to the half-maximum of the peak in the elimination phase) to the whole data area because a peak function could not cover the whole data. It indicated that the signals required two or more functions to fit the whole data. Furthermore, the tail part should be applied with the same function in the peak area because the circulation system could not be changed while the signal was measured. Therefore, a high possibility existed that the whole data curve of the single peak would be fitted by linear bi-summations of the peak functions. Fig. 2B presents the fitting results of the five peak functions. All linear bi-summation functions fitted as shown in Fig. 2C. To compare the results, we calculated the R-squared values as shown in Table 1. We excluded the result of the bi-Gaussian function because it was too deviated to fit on the curve.

The clinical data provided the best fit to the EMG function on both the peak and whole data area (Table 1). This result demonstrated that the multi-EMG function was required for fitting on multiple peak curves. The same number of EMG functions was needed; the peak and tail parts required only one.

### Multi-EMG modelling for ICG curve

The bi-EMG function was well suited to expressing the single peak data with the time profile. We verified whether the multi-EMG function would fit on various forms of the signals. We collected several quantitative ICG-VA signals, which had a couple of peaks and a tail intact and varied depending on the status of the blood flow. Furthermore, we evaluated the multi-EMG model as the universal fitting function for the cerebrovascular hemodynamic curve (Fig. 3). The signals were fitted to ter- or quad-EMG, specifically three or four respective linear summations of EMG. The R-squared values revealed that the multi-EMG functions fitted well with the various forms of data (average value, 0.9832).

In addition, we identified the center of peaks at the peak area of the fitting curve and measured data using the peakfinding algorithm with the first derivative method. The first derivative values of the given function were equal to zero at the local peak point. The low degree Savitzky-Golay filter, one of the smoothing filters, was applied to the measured data while finding the center of peaks. The location of the center of peaks at the time axis in the fitting curve (grey numbers) was similar to that of the measured data (black numbers). Therefore, we could distinguish the center of peaks in the clinical data from one in the multi-EMG function.

The partial data were sometimes missed, which may have been occurred when the intensity of the signal exceeded the measurement system limit. The fitting model restored the parameters of the particular missing data through the data trend. Fig. 4 shows that the multi-EMG model recovered the missing data and predicted center of peaks that could not be identified on the original measured data. To identify the predicted center of peaks, we used the same method with the fitting functions of non-saturated data.

## DISCUSSION

In the evaluation, the numerical model of the quantitative ICG-VA provided accurate information on cerebrovascular hemodynamics and enabled a comparison of the blood flow in various conditions by using the primary parameters. The multi-EMG function was the best numerical function for reflecting cerebrovascular hemodynamics, which was proved using a nonlinear regression fitting algorithm. The number of EMG functions in multi-EMG was determined by the number of peaks in the clinical data. This is because the peak and tail were required for each component of the EMG function. For each peak, one component of the EMG function was dominant.

Fig. 5A shows the two component EMG functions composed of a single peak fitting function. The peak area of the fitting curve is predominantly influenced by the growth part of the first EMG. In terms of the logarithmic scale of the curve (Fig. 5B), the decay part of the peak area aligns with the exponential line in the first EMG. As shown in these figures, the first EMG component is dominant in the peak area; the second strongly influences the tail part. Therefore, we considered a relevant single EMG to analyse the properties of each peak.

An objective of curve fitting with nonlinear regression is to set a numerical model describing a complex biological system using equations and numerical parameters. The numerical models enable interpretation of the measured data as understandable systems and estimation of the parameters or trends of curves. Moreover, if data are missing from the “valid part”—which may be caused by the determination constraint in the measurement system—the model not only predicts the parameter, but it also restores the missing data to some extent [16].

Thus, to evaluate the model’s effectiveness in data restoration, we maintained the top of the main peak from 60% to 100% in the real measured data based on the peak height (Fig. 6A and C). In the simulation, the intensity of the measurement signals exceeded the observable maximum of the ICG-VA system, as shown in Fig. 6. We reconstructed the missing parts at each truncated main peak (Fig. 6B and D) by applying the multi-EMG model. In addition, we estimated the primary parameters (center of peaks, peak height, steepness of growth (w), and the decay constant (

Fig. 6E presents the relative difference based on the parameters of the original fitting function. The accuracy of the reconstruction of all parameters is usually lower when the larger missing parts exist; however, the accuracy has a <30% margin of error if the peak remains greater than 70%. The decay constant has a relatively lower accuracy than the growth steepness because the larger missing parts influence the second EMG located in the tail part. However, the decay constant values are reliable if there are small missing areas, making it a possible evaluation value for the properties of the elimination phase.

There are other methods for assessment of hemodynamics using ICG-VA. Flow 800 (Carl Zess, Oberkochen, Germany) is a well-known software and there are several reports on its use in cerebrovascular surgery [6,32]. It is convenient for surgeons, without requiring any calculation by users. However, because the software can only calculate the delay and the slope (average intensity, AI/second) at a specific point of the intensity graph within the range arbitrarily determined by the clinician or researcher, there is a limit to obtain accurate quantitative values [7]. Our method provides the coefficient of fitting function based on nonlinear regression analysis, and it offers more objective evaluation indicators. Moreover, the commercial software is an integrated part of the operative microscope, compatible models were launched on 2011 and of considerable price (approximately 30 million Korean Won). Besides, it is impossible for analysis of data acquired in older model of microscope.

Considering the pros and cons of the commercial software and our methods, it might be said : commercial software is useful for intraoperative assessment for hemodynamics on site, whereas our methods for analysis for data in batch mode, especially for archived data to elucidate hemodynamic features of series of cases in more objective way with better accuracy.

The above results showed that the multi-EMG model reconstructed missing peak parameters with reliable prediction in the smaller missing area. The predicted values of the center of peaks in Fig. 6 are acceptable for the quantitative analysis of ICG-VA.

## CONCLUSION

Establishing a numerical model for analysing a biological system has several advantages [25,26]. First, a mathematical model enables visualisation and intuitive comprehension of biological phenomena [30]. Biological phenomena usually occur as a result of combined reactions to a complex system. The numerical model developed in this study simplifies a given phenomenon into several equations and parameters. Through an analysis of these equations and parameters, the complex system can then be understood. Second, the trend and response of a phenomenon can be predicted as presented by the restoration of data in our series. Furthermore, the tempolabile quantitative analysis presents a time profile of vascular blood flow and enables a more accurate evaluation.

Based on the results of this study, the multi-EMG model can be used as a mathematical model for cerebrovascular hemodynamics to understand vascular blood flow, which is quantitatively measured by ICG-VA during cerebrovascular surgery. The proposed multi-EMG model fitted well in various conditions of vascular blood flow by nonlinear regression and reflected the parameters (such as center of peaks) of clinical data. Furthermore, we proved that the multi-EMG model reconstructed a portion of missing data and predicted the parameters of these missing parts. Our numerical model can thus enhance the insight of cerebrovascular hemodynamics in various clinical conditions.

## Notes

**Conflicts of interest**

Hee-Jin Yang has been editorial board of JKNS since May 2020. He was not involved in the review process of this original article. No potential conflict of interest relevant to this article was reported.

**Informed consent**

This type of study does not require informed consent.

**Author contributions**

Conceptualization : HJY, HC; Data curation : HJY, HC, PSS; Formal analysis : HJY, HC, PSS; Funding acquisition : HJY, JHS; Methodology : HJY, HC, PSS; Project administration : HJY; Visualization : HJY, HC; Writing - original draft : HC, YJS, SBP, PSS, JHS, HJY; Writing - review & editing : HC, YJS, SBP, PSS, JHS, HJY

**Data sharing**

The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.

**Preprint**

None

## Acknowledgements

This work was partly supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2021R1F1A1056527). This work was partly supported by 2022 University of Seoul Research Grant. This work was partly supported by a clinical research fund of Seoul National University Boramae Hospital.