ADVANCED NON-INVASIVE PHOTONIC METHODS FOR FUNCTIONAL MONITORING OF HAEMODYNAMICS AND VASOMOTOR REGULATION IN HEALTH AND DISEASES

EDITED BY : Alexey Goltsov, Viktor V. Sidorov, Sergei G. Sokolovski and Edik Rafailov PUBLISHED IN : Frontiers in Physiology

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-760-7 DOI 10.3389/978-2-88963-760-7

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# ADVANCED NON-INVASIVE PHOTONIC METHODS FOR FUNCTIONAL MONITORING OF HAEMODYNAMICS AND VASOMOTOR REGULATION IN HEALTH AND DISEASES

Topic Editors:

Alexey Goltsov, Russian Technological University - MIREA, Russia Viktor V. Sidorov, SPE "LAZMA" Ltd., Russia Sergei G. Sokolovski, Aston University, United Kingdom Edik Rafailov, Aston University, United Kingdom

Citation: Goltsov, A., Sidorov, V. V., Sokolovski, S. G., Rafailov, E., eds. (2020). Advanced non-invasive photonic methods for functional monitoring of haemodynamics and vasomotor regulation in health and diseases. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-760-7

# Table of Contents

*05 Editorial: Advanced Non-invasive Photonic Methods for Functional Monitoring of Haemodynamics and Vasomotor Regulation in Health and Diseases*

Alexey Goltsov, Viktor V. Sidorov, Sergei G. Sokolovski and Edik U. Rafailov

*07* Matrix metalloproteinase*-9 Gene-1562C>T Gene Polymorphism and Coronary Artery Disease in the Chinese Han Population: A Meta-Analysis of 5468 Subjects*

Yan-Yan Li, Xin-Xing Yang, Yan-Hong Zhou, Ge Gong, Hong-Yu Geng, Hyun J. Kim, Chuan-Wei Zhou, Yun Qian, Xiang-Ming Wang and Jun Wu


Yurie Obata, Qi J. Ong, J. T. Magruder, Helen Grichkevitch, Dan E. Berkowitz, Daniel Nyhan, Jochen Steppan and Viachaslau Barodka

*52 Bifurcation in Blood Oscillatory Rhythms for Patients with Ischemic Stroke: A Small Scale Clinical Trial Using Laser Doppler Flowmetry and Computational Modeling of Vasomotion*

Alexey Goltsov, Anastasia V. Anisimova, Maria Zakharkina, Alexander I. Krupatkin, Viktor V. Sidorov, Sergei G. Sokolovski and Edik Rafailov

*63 Physical and Chemical Processes and the Morphofunctional Characteristics of Human Erythrocytes in Hyperglycaemia* Victor V. Revin, Natalia A. Klenova, Natalia V. Gromova, Igor P. Grunyushkin, Ilia N. Solomadin, Alexander Y. Tychkov, Anastasia A. Pestryakova,

Anna V. Sadykhova, Elvira S. Revina, Ksenia V. Prosnikova, Jean-Christophe Bourdon and Nikolai Zhelev


Wei-tong Guo, Hongyin Ma, Jia Liu, Zhen-Ni Guo and Yi Yang

#### *88 Multimodal Optical Diagnostics of the Microhaemodynamics in Upper and Lower Limbs*

Angelina I. Zherebtsova, Viktor V. Dremin, Irina N. Makovik, Evgeny A. Zherebtsov, Andrey V. Dunaev, Alexey Goltsov, Sergei G. Sokolovski and Edik U. Rafailov

#### *103 Oligonucleotide Microarrays Identified Potential Regulatory Genes Related to Early Outward Arterial Remodeling Induced by Tissue Plasminogen Activator*

Olga Plekhanova, Yelena Parfyonova, Irina Beloglazova, Bradford C. Berk and Vsevolod Tkachuk

# Editorial: Advanced Non-invasive Photonic Methods for Functional Monitoring of Haemodynamics and Vasomotor Regulation in Health and Diseases

Alexey Goltsov 1,2 \*, Viktor V. Sidorov <sup>3</sup> , Sergei G. Sokolovski <sup>4</sup> and Edik U. Rafailov <sup>4</sup>

*1 Institute of Cybernetics, Russian Technological University, Moscow, Russia, <sup>2</sup> M&S Decisions LLD, Moscow, Russia, <sup>3</sup> SPE "LAZMA" Ltd., Moscow, Russia, <sup>4</sup> Optoelectronics and Biomedical Photonics Group, Aston Institute of Photonic Technologies, Aston University, Birmingham, United Kingdom*

Keywords: non-invasive multimodal photonic diagnostics, laser Doppler flowmetry, NIRS, optoacoustic technique, ischemic stroke, hemodynamics regulation

**Editorial on the Research Topic**

#### **Advanced Non-invasive Photonic Methods for Functional Monitoring of Haemodynamics and Vasomotor Regulation in Health and Diseases**

Edited and reviewed by: *Claudia Penna, University of Turin, Italy*

\*Correspondence: *Alexey Goltsov alexey.goltsov@gmail.com*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *08 February 2020* Accepted: *20 March 2020* Published: *15 April 2020*

#### Citation:

*Goltsov A, Sidorov VV, Sokolovski SG and Rafailov EU (2020) Editorial: Advanced Non-invasive Photonic Methods for Functional Monitoring of Haemodynamics and Vasomotor Regulation in Health and Diseases. Front. Physiol. 11:325. doi: 10.3389/fphys.2020.00325* In the present Research Topic we collected papers devoted to the different aspects of the noninvasive photonic techniques capable monitoring real-time microvascular hemodynamics in health and pathology. Compact laser instruments for near-infrared spectroscopy (NIRS), Laser Doppler Flowmetry (LDF), diffuse correlation spectroscopy (DCS), NIRS flow oximetry, and others are routinely and widely used for concurrent non-invasive measurement of different parameters of hemocirculation such as blood perfusion, mean arterial pressure, oxy- and deoxyhemoglobin concentrations, and others. Further development and introduction of these approaches and systems to clinical practice will lead to enrichment and step-by-step substitution of laboratory tests with enabling bed-side diagnostics and wireless distantly controlled patient condition monitoring.

The biophotonic approaches became a powerful tool to study dynamic processes such as vascular tone autoregulation and functional hyperaemia which are controlled by vessel vasculature rhythms derived from different regulatory systems: nervous sympathetic and parasympathetic, hormonal and mechanical muscle regulation. A promising area for research is establishing an autoregulatory link between blood flow and level of the metabolic activity where combined LDF and fluorescent approaches can bring significant advantages in the diagnostics. Spectral analysis of different optical data arrays can provide valuable information on neuronal-vasculature and muscular-vasculature coupling autoregulation mechanism under intense stress and/or physical activity in health and at various pathologies. Due to the morphological and physiological differences in the microcirculation bed of the human skin, a non-invasive detailed research of blood microcirculation can be organized using a new class of wearable wireless laser flowmeters adopted to the range of the research areas.

This Research Topic brings together the cutting-edge applications of non-invasive biophotonic approaches and techniques allowing to translate personalized data array analysis to distinguished diagnosis.

In the research paper, Goltsov et al. applied spectral analysis of the LDF signals for the investigation of cerebrovascular hemodynamics in patients with post-acute ischemic stroke. The observed asymmetry in LDF signals measured on the stroke-affected and unaffected hemispheres of the patients was analyzed by computational modeling of vasomotion of arteriole diameter.

Zherebtsova et al. discussed the development and application of the optical non-invasive diagnostic approach for the detection and evaluation of the microcirculatory severity and metabolic disorders in patients with rheumatic diseases and diabetes mellitus. The proposed methods include the joint use of LDF, absorption spectroscopy and fluorescence spectroscopy in combination with functional tests in the field of rheumatology and endocrinology.

Semyachkina-Glushkovskaya et al. investigated brain hemorrhages in newborn rats using an integrative approach including magnetic resonance imaging, monitoring of cerebral blood flow using Doppler coherent tomography, oximetry, coherent-domain optical technologies for visualization of the cerebral blood flow, investigations of the deformability of red blood cells, and morphological analysis of brain tissues. Combination of these techniques in the monitoring of the preand post-hemorrhage periods allowed authors to investigate a time-dependent stress inducing changes in neonatal brain areas which are involved in higher cognitive functions in progress of the pathology.

In the method article, Caicedo et al. developed a new computational method for the decomposition of NIRS signals into the partial linear contributions of different physiological signals. Authors combined Oblique Subspace Projections in order to define adequate signal subspaces and the linear regression model of NIRS signals to calculate regression parameters and estimate the contribution of various input variables. Authors evaluate the performance of the developed methods by case studies in the field of cerebral hemodynamics monitoring using NIRS and concomitant measurements of MABP, EtCO2, HR, SaO2, and ECMO signals.

In mini-survey, Esenaliev reviewed optoacoustic technique which is a novel medical diagnostic platform for non-invasive measurement of physiologic variables, functional imaging, and hemodynamics monitoring. Authors discussed the application of optoacoustic methods to the measurement of important physiologic parameters including temperature, thermal coagulation, freezing, concentration of molecular dyes and nanoparticles in organs, oxygenation, and hemoglobin concentration with high temporal and spatial resolution.

Guo et al. analyse the daily rhythm of dynamic cerebral autoregulation by simultaneous measurement of the arterial blood pressure at the brachial artery and cerebral blood flow velocity over the course of a day in healthy adults.

Obata et al. used simultaneously recorded the ECG and plethysmograph in the investigation of the effect of body position and exercise on pulse arrival time taken for the pulse wave to reach the vasculature of various peripheral tissue beds. Authors analyzed vascular autoregulation mechanisms that maintain the optimal arrival time irrespective of its distance from the heart.

Several papers in the Research Topic are devoted to the biophysical investigation of vascular regulation at the molecular, cellular, and tissue levels. Revin et al. presented results of the detailed biophysical investigation on the hyperglycaemia effect on morphological changes in erythrocytes. Authors combined Raman spectroscopy in the observation of conformational changes in hemoglobin with laser interference microscopy to determine structural transition and a composition change in the erythrocyte membranes at graduated hyperglycaemia. Ivanov et al. investigated the possibility of the application of microwave radiometry method to non-invasive monitoring of the physicochemical properties of circulating proteins. In brief research report article, authors provided results on the application of commercially available microwave radiometry equipment to the investigation of denaturation kinetics of albumin in solution. Plekhanova et al. investigated the potential role of plasminogen activators (PA) in the promotion of outward arterial remodeling after injury in rat common carotids. Authors studied a gene expression profile change after injury and local treatment with the recombinant PA in the vessel wall related to vascular tone. Li et al. discussed the prospective researches of the plasma matrix metalloproteinase-9 (MMP-9) contents that can serve as a predictive indicator for cardiovascular disease mortality in the Chinese Han population.

The results presented in this Research Topic by different research groups show a wide spectrum of the current directions for the development and application of biphotonic techniques in medicine. As seen, the mainstream research in biophotonics focuses on the monitoring and processing of concurrently collected biomedical signals in patients to enhance current clinical diagnosis and develop novel diagnostic techniques and equipment.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

**Conflict of Interest:** AG is employed by company M&S Decisions LLD. VS is employed by the company SPE "LAZMA" Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Goltsov, Sidorov, Sokolovski and Rafailov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Matrix metalloproteinase-9 Gene-1562C>T Gene Polymorphism and Coronary Artery Disease in the Chinese Han Population: A Meta-Analysis of 5468 Subjects

Yan-Yan Li <sup>1</sup> \*, Xin-Xing Yang<sup>1</sup> , Yan-Hong Zhou<sup>1</sup> , Ge Gong<sup>1</sup> , Hong-Yu Geng<sup>1</sup> , Hyun J. Kim<sup>2</sup> , Chuan-Wei Zhou<sup>1</sup> , Yun Qian<sup>1</sup> , Xiang-Ming Wang<sup>1</sup> and Jun Wu<sup>1</sup>

*<sup>1</sup> Department of Geriatrics, First Affiliated Hospital of Nanjing Medical University, Nanjing, China, <sup>2</sup> Department of Radiation Oncology, University of Pennsylvania, Philadelphia, PA, USA*

Background: Multiple studies indicate that the *matrix metalloproteinase-9 (MMP-9)*-1562C>T gene polymorphism may be associated with an increased risk of coronary artery disease (CAD) in the Chinese Han population. However, a clear consensus has yet to be established.

#### Edited by:

*Alexey Goltsov, Abertay University, UK*

## Reviewed by:

*Karen Yvonne Stokes, Louisiana State University Health Sciences Center, USA Wang Min, Yale University, USA*

> \*Correspondence: *Yan-Yan Li lyynjmu123@126.com*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *30 December 2015* Accepted: *23 May 2016* Published: *09 June 2016*

#### Citation:

*Li Y-Y, Yang X-X, Zhou Y-H, Gong G, Geng H-Y, Kim HJ, Zhou C-W, Qian Y, Wang X-M and Wu J (2016) Matrix metalloproteinase-9 Gene-1562C*>*T Gene Polymorphism and Coronary Artery Disease in the Chinese Han Population: A Meta-Analysis of 5468 Subjects. Front. Physiol. 7:212. doi: 10.3389/fphys.2016.00212* Objective and methods: A meta-analysis of 5468 subjects from 10 separate studies was performed to explore the possible relationship between the *MMP-9*-1562C>T gene polymorphism and CAD within the Chinese Han population. Pooled odds ratio (ORs) for the association and the corresponding 95% confidence intervals (CIs) were evaluated by a random or fixed-effect model.

Results: Our analysis confirms the association between the *MMP-9*-1562C>T gene polymorphism and an increased risk of CAD within the Chinese Han population under allelic (OR: 1.60, 95% CI: 1.25–2.04, *P* = 0.0002), recessive (OR: 3.05, 95% CI: 1.67–5.56, P = 0.0003), dominant (OR: 2.23, 95% CI: 1.49–3.35, *P* = 0.0001), homozygous (OR: 3.41, 95% CI: 1.87–6.23, *P* < 0.0001), heterozygous (OR: 2.03, 95% CI: 1.40–2.93, *P* = 0.0002), and additive genetic models (OR: 1.78, 95% CI: 1.33–2.39, *P* < 0.0001).

Conclusions: In the Chinese Han population, the *MMP-9*-1562C>T gene polymorphism is correlated with an increased risk of CAD. Therefore, Han Chinese carriers of the -1562T allele may be at an increased risk of CAD.

Keywords: matrix metalloproteinase-9, -1562C>T, polymorphism, coronary artery disease, Chinese

#### INTRODUCTION

Coronary artery disease (CAD) is a chronic condition with both hereditary and environmental factors. The increasing prevalence of unhealthy lifestyles in China (i.e., extended periods of sitting, a more sedentary lifestyle, a high-fat diet), the general aging of the Chinese population, and the increasing rates of hypertension and diabetes have damaged the quality of life of many Chinese

**7**

people and sharply increased CAD-associated morbidity and mortality (Li, 2007). However, recent advancements in molecular biology have allowed researchers to better elucidate the pathogenic mechanism of CAD.

Matrix metalloproteinases (MMPs), a large family of zincdependent proteolytic enzymes that have special functions in extracellular matrix (ECM) degradation have been a particular rich area of research. MMPs play a role in multiple pathophysiological processes such as ECM degradation, inflammation response, tumor metastasis, and atherosclerosis (Castro and Tanus-Santos, 2013). Enhanced MMP expression and activity also play a pivotal role in early arteriosclerosis, plaque rupture, myocardial infarction, and heart failure (Tanner et al., 2011). MMP-9 (also known as gelatinase B) is one the most extensively researched MMPs and exhibits enhanced expression in atherosclerosis injury sites. It has been also found that there were many gene sequence variations in MMP-9 gene, of which the locus rs3918242 (-1562C>T) in the promoter region was the most reported gene locus.

MMP-9 gene, located in 20q11.2–13.1, spans 7.7 kb and contains 13 exons. MMP-9-1562C>T mutation occurs in a substitution of cytosine (C) for thymine (T). This mutation may increase MMP-9 gene expression and risk of CAD development by reducing the binding affinity with transcription inhibition proteins.

Although many studies on the association between MMP-9- 1562C>T gene polymorphism and CAD have been performed in China, the composite of these studies fail to provide a consensus. In 2005, Tang et al found that MMP-9-1562C>T gene polymorphism was significantly associated with CAD in a Zhejiang population and T allele increased the CAD risk (Tang et al., 2005). In 2007, Chen et al also reported a similar result in another Hunan population (Chen et al., 2007). In contrast, Wu et al failed to find an association between MMP-9-1562C>T gene polymorphism and CAD in a Beijing population in 2009 (Wu et al., 2009) and Zhi et al observed no significant effects for MMP-9-1562C>T gene polymorphism on CAD risk in a Jiangsu population (Zhi et al., 2010).

We performed the current meta-analysis with the hopes to provide a valuable conclusion on the association between MMP-9-1562C>T gene polymorphism and CAD in the Chinese Han population.

# MATERIALS AND METHODS

#### Publication Search and Inclusion Criteria

The Web of Science, PubMed, Embase, the China Biological Medicine Database, and the China National Knowledge Infrastructure electronic databases were searched using the terms "matrix metalloproteinase-9," "coronary artery disease," or "coronary heart disease," and "polymorphism" in our initial search. Retrieved studies were published between 2005 and 2010 with the last study updated on March 26, 2016.

To meet our inclusion criteria, studies had to (a) evaluate the association between MMP-9-1562C>T gene polymorphism and CAD in the Chinese Han population. (b) diagnose CAD according to the clinical symptoms combined with examination results (i.e., coronary arteriography, electrocardiogram, treadmill exercise test, echocardiography, myocardial perfusion imaging by Emission Computed Tomography, etc.) with a minimal stenosis rate of the major coronary artery diameter of more than 50% (c) be a case-control or cohort study published in an official journals or as a postgraduate dissertation.

# Data Extraction

Data was extracted according to a standardized protocol. Two investigators searched for duplicates while a third served as an arbiter to resolve possible disagreements. Duplicate papers, those that violated the inclusion criteria, or those that provided deficient data were removed. Identical data sets used in different studies by the same authors were used once. Abstracted data consisted of the following items: the first author's name, publication year, region, number of genotypes, genotyping method, study design, age, gender, and total number of cases and controls.

# Statistical Analysis

The odds ratio (OR) and its corresponding to 95% confidence interval (CI) were used to measure the association between MMP-9-1562C>T gene polymorphism and CAD. The Chisquare-based Q-test was adopted to measure the effects betweenstudies heterogeneity (P < 0.05 level) (Cochran, 1968).Variation due to heterogeneity was assessed by the inconsistency index I 2 . If heterogeneity were present in the study, the randomeffects model would be used to assess the combined OR (the DerSimonian and Laird method) (DerSimonian and Laird, 1986). Otherwise, the fixed-effects model would be used (the Mantel-Haenszel method) (Mantel and Haenszel, 1959). The pooled OR was determined by Z-test and significance was set at P < 0.05.

The Fisher's exact test was used to evaluate the Hardy-Weinberg equilibrium (HWE) (P < 0.05). Potential publication bias was assessed by Egger′ s linear regression test on the natural log scale of the OR to detect funnel plot asymmetry (P < 0.05 level) (Egger et al., 1997). The statistical analysis was performed by STATA 11.0 software (StataCorp, College Station, TX).

# RESULTS

## Studies and Populations

Ten of the nineteen retrieved papers fit the inclusion criteria. Of the nine excluded studies, one was a duplicate, four were reviews, and another four were irrelevant to our interests. No study was excluded for deviation from HWE. The ten studies compiled the data from 3168CAD patients and 2300 controls (**Table 1**, Presentation 1 in Supplementary Material; Tang et al., 2005; Meng et al., 2006; Chen et al., 2007; Wang et al., 2007; Wu et al., 2009; Gao and Wang, 2010; Ma et al., 2010; Yong and Shi, 2010; Zhang et al., 2010; Zhi et al., 2010) and represented seven provinces (Shanxi, Xinjiang, Jiangsu, Zhejiang, Hunan, Tianjin, and Beijing). All subjects were of Han ethnicity. However, there were substantial differences in the size of the total patient population between these studies. The Wu N study alone contributed 1356 CAD patients and 689 controls for this


TABLE 1 | Characteristics of the investigated studies of the association between the matrix metalloproteinase-9 gene -1562C>T polymorphism and coronary artery disease in the Chinese population.

*CAD, coronary artery disease.*

*Polymerase chain reaction-restriction fragment length polymorphism genotyping method and Case-control study design were adopted in the above studies.*

*Compared with control group,* \**P* < *0.05.*

TABLE 2 | Summary of meta-analysis of association of matrix metalloproteinase-9 gene -1562C>T polymorphism and coronary artery disease in the Chinese population.


\**P* < *0.05. CAD, coronary artery disease; CI, confidence interval; OR, odds ratio; CAD size, the total number of CAD cases; control size, the total number of control group; homo genetic model, homozygous genetic model; hetero genetic model, heterozygous genetic model; CT1, CT sample size of CAD group; CC0, CC sample size of control group; T1, Total sample size of CAD group.*

meta-analysis. Differences in patient number may be factor in the lack of consensus on this topic.

## Combined Analyses

There was a significant association between MMP-9-1562C>T gene polymorphism and CAD in the Chinese Han population under allelic (OR: 1.60, 95% CI: 1.25–2.04, P = 0.0002), recessive (OR: 3.05, 95% CI: 1.67–5.56, P = 0.0003), dominant (OR: 2.23, 95% CI: 1.49–3.35, P = 0.0001), homozygous (OR: 3.41, 95% CI: 1.87–6.23, P < 0.0001), heterozygous (OR: 2.03, 95% CI: 1.40–2.93, P = 0.0002), and additive genetic models (OR: 1.78, 95% CI: 1.33–2.39, P < 0.0001). (**Table 2**, **Figures 1**–**6**).

There was also significant heterogeneity under the allelic (P = 0.002, I <sup>2</sup> = 65.1%), dominant (P < 0.00001, I <sup>2</sup> = 83.5%), heterozygous (P < 0.00001, I <sup>2</sup> = 80.3%) and additive genetic models (P < 0.0001, I <sup>2</sup> = 75.9%). Subsequent meta-regressions explored the source of this heterogeneity under their respective genetic models.

Under the allelic genetic model, CC sample size of CAD group (CC1, P = 0.006), CT sample size of CAD group (CT1, P = 0.005), TT sample size of CAD group (TT1, P = 0.010), CC sample size of control group (CC0, P = 0.023), total sample size of control group (T0, P = 0.034), and weight (%) (P = 0.030)were possible confounding factors that could partially explain heterogeneity between studies.

Under the additive genetic model, CC1 (P = 0.013), CT1 (P = 0.012), and TT1 (P = 0.016) were the confounding factors. In the allelic and additive genetic models, CT1 plays a central role


FIGURE 1 | Forest plot of coronary artery disease associated with MMP-9-1562C>T gene polymorphism under an allelic genetic model stratified by CT1 (distribution of T allelic frequency of MMP-9 gene).

CC+CT).

in explaining the source of heterogeneity. According to CT1, the whole population was separated into two subgroups. The studies with CT1 < 30 were grouped into subgroup 1 and the residual studies with CT1 ≥ 60 belonged to subgroup 2.

Stratified by CT1, a subgroup analysis under allelic and additive genetic models found significant association between MMP-9-1562C>T gene polymorphism and CAD in both subgroups (allelic genetic model: subgroup 1: OR: 1.74, 95% CI: 1.28–2.36, P = 0.0004; subgroup 2: OR: 1.52, 95% CI: 1.10–2.11, P = 0.01) (additive genetic model: subgroup 1: OR: 1.89, 95% CI: 1.34–2.68, P = 0.0003; subgroup 2: OR: 1.72, 95% CI: 1.14–2.57, P = 0.009). No significant


FIGURE 3 | Forest plot of coronary artery disease associated with MMP-9-1562C>T gene polymorphism under a dominant genetic model stratified by CC0 (CT+TT vs. CC).

heterogeneity was fond in subgroup 1 (allelic genetic model: Pheterogeneity = 0.38, I <sup>2</sup> = 4.0%; additive genetic model: Pheterogeneity = 0.26, I <sup>2</sup> = 23.7%), but significant heterogeneity was detected in subgroup 2 (allelic genetic model: Pheterogeneity = 0.001, I <sup>2</sup> = 78.3%; additive genetic model: Pheterogeneity < 0.0001, I <sup>2</sup> = 85.8%). This identifies CT1 as the primary confounding factor under the allelic and additive genetic models (**Tables 2**–**4**; **Figures 1**, **6**).

Under the dominant genetic model, CC1 (P = 0.015), CT1 (P = 0.009), TT1 (P = 0.009), CC sample size of control group (CC0, P = 0.001), CT sample size of control group (CT0, P = 0.001), TT sample size of control group (TT0, P = 0.012), and


FIGURE 5 | Forest plot of coronary artery disease associated with MMP-9-1562C>T gene polymorphism under a heterozygous genetic model stratified by T1 (CT vs. CC).


FIGURE 6 | Forest plot of coronary artery disease associated with MMP-9-1562C>T gene polymorphism under an additive genetic model stratified by CT1 (Total T vs. Total C).


TABLE 3 | The meta-regression results among 10 studies in the Chinese population under an allelic genetic model for matrix metalloproteinase-9 gene -1562C>T gene polymorphism.

\**P* < *0.05. Coefficient: regression coefficient. The regression coefficients are the estimated increase in the lnOR per unit increase in the covariates. Cons, constant item.*

weight (P = 0.001) could partly explain the heterogeneity source. In a subgroup analysis stratified by CC0, significant increased CAD risk was observed in both subgroups (subgroup 1: OR: 1.73, 95% CI: 1.13–2.65, P = 0.01; subgroup 2: OR: 3.40, 95% CI: 1.41–8.21, P = 0.006). Although heterogeneity was still detected in both subgroups, it was significantly reduced in subgroup 2 (Pheterogeneity = 0.0003, I <sup>2</sup> = 80.9%), suggesting that CC0 was the main confounding factor (**Tables 2, 5**; **Figure 3**).

Under the heterozygous genetic model, CC1 (P = 0.006), total sample size of CAD group (T1, P = 0.011), and CT0 (P = 0.004) could partly explain the heterogeneity. In a subgroup analysis stratified by T1, significant increase in CAD risk was detected in both subgroups (subgroup 1: OR: 3.19, 95% CI: 1.95–5.21, P < 0.00001; subgroup 2: OR: 1.41, 95% CI: 1.00–2.01, P = 0.05). No significant heterogeneity existed in subgroup 1 any longer (Pheterogeneity = 0.12, I <sup>2</sup> = 45.5%), but significant heterogeneity was still observed in subgroup 2 (Pheterogeneity < 0.0001, I <sup>2</sup> = 85.8%; **Tables 2, 6**; **Figure 5**).

### Sensitivity Analysis

In the current meta-analysis, the sensitivity analysis was performed. After the Wu et al. was removed from the current meta-analysis, the association between MMP-9-1562C>T gene polymorphism and CAD was further strengthened under the allelic genetic model (OR: 1.72, 95% CI: 1.35–2.19, P = 1.3 × 10−<sup>5</sup> , Pheterogeneity = 0.04, I <sup>2</sup> = 50.0%; Wu et al., 2009). Removal of other studies respectively from the current studies did not change the results from the original analysis. Hence, the Wu et al study should be the high sensitivity study in the current meta-analysis.

#### Bias Diagnostics

The publication bias of the studies was evaluated by funnel plot and Egger's test. There was no visual publication bias in the funnel plot (**Figure 7**). No statistically significant difference was detected in the Egger's test, implying no publication bias existed in the current meta-analysis under the recessive genetic model (T = −0.73, P = 0.487).

# DISCUSSION

In the present meta-analysis, we found a significant association between MMP-9-1562C>T gene polymorphism and CAD in the Chinese Han population under allelic (OR: 1.60), recessive (OR: 3.05), dominant (OR: 2.23), homozygous (OR: 3.41), heterozygous (OR: 2.03), and additive genetic models (OR: 1.78). Hence, it has been concluded that in Chinese Han population, the MMP-9-1562C>T gene polymorphism may be associated with the increased CAD susceptibility among Han Chinese.

What contributed to the recent controversy over the association between MMP-9-1562C>T gene polymorphism and CAD? The meta-regression used to reveal the source of the heterogeneity detected under allelic, dominant, heterozygous and additive genetic models (Pheterogeneity < 0.05) suggests that patient number may have been the confounding factor. In the heterogeneity source analysis, CT1 was possibly indicated to be the main heterogeneity source under allelic and additive genetic models. Although the subgroup analysis stratified by CT1 showed a significantly increased risk of CAD in both subgroups, only one subgroup exhibited heterogeneity. Hence, CT1 was the main confounding factor contributing to the heterogeneity under the allelic and additive models. Similarly, T1 and CC0 were major confounding factors under heterozygous and dominant models, respectively. CT1, T1, CC0 may be better matched between the CAD and control groups under these genetic models.

Sensitivity analysis showed that our meta-analysis was most sensitive to the Wu et al study. Although the pooled analysis result without Wu et al study was different from the original result (Wu et al., 2009), they are still consistent in associating the presence of the gene polymorphism with an increased risk of CAD.

MMPs belong to a neutral protease family that contains zinc ions. MMP-9 is the leading MMP expressed and secreted by the vascular cell walls. It is also secreted by monocytes, neutrophils, vascular smooth muscle cells (VSMCs), and endothelial cells. The relative molecular weight of MMP-9 is 92KD. In its active form its molecular weight is 84KD. MMP-9 can degrade extensive ECM substrates, including Type IV collagen, which plays a key role in the revascularization, inflammation response, and atherosclerosis progression. Research has shown an increased MMP-9 expression level in the atherosclerotic arteries of human and animals compared to normal arteries. The MMP-9 degradation activity was most located in the shoulder regions of plaque, the lipid core margin, and micro-vessels formation regions. This suggests that MMP-9 may be associated with the coronary artery TABLE 4 | The meta-regression results among 10 studies in the Chinese population under an additive genetic model for matrix metalloproteinase-9 gene -1562C>T gene polymorphism.


\**P* < *0.05. Coefficient: regression coefficient. The regression coefficients are the estimated increase in the lnOR per unit increase in the covariates. Cons, constant item.*

TABLE 5 | The meta-regression results among 10 studies in the Chinese population under a dominant genetic model for matrix metalloproteinase-9 gene -1562C>T gene polymorphism.


\**P* <*0.05. Coefficient: regression coefficient. The regression coefficients are the estimated increase in the lnOR per unit increase in the covariates. Cons, constant item.*

TABLE 6 | The meta-regression results among 10 studies in the Chinese population under a heterozygous genetic model for matrix metalloproteinase-9 gene -1562C>T gene polymorphism.


\**P* < *0.05. Coefficient: regression coefficient. The regression coefficients are the estimated increase in the lnOR per unit increase in the covariates. Cons, constant item.*

plaque stability and myocardial infarction (Speidl et al., 2011).

The animal experiments have discovered that the atherosclerosis lesions and VSMCs intima migration in the MMP-9 gene knock-out mice were remarkably decreased than that in the wild-type mice (Ye, 2006). The clinical researches have confirmed that the high MMP-9 expression level was correlated with the premature CAD, unstability of the coronary atherosclerosis plaque, the in-stent restenosis and arterial aneurysm formation (Jones et al., 2006). The prospective researches have shown that plasma MMP-9 contents can serve as the prediction indicator for the cardiovascular diseases mortality (Blankenberg et al., 2003).

MM9 is regulated primarily at the transcriptional level. The MMP-9-1562C>T gene polymorphism is located in crucial regulatory elements, including the 9 bp sequence GCGCAC/TGCC (−1567→−1559), a potential binding site for transcription inhibition proteins (Zhang et al., 1999). In 2002, Cho et al reported a change in bond zone structure and a weakened binding capacity between DNA and transcription inhibition protein, when the MMP-9-1562C allele was replaced by -1562T allele (Cho et al., 2002). This generates a high and low activity promoter genotype (CT/TT and CC, respectively) with increased transcription with the high activity promoter genotypes.

Increased MMP-9 expression may contribute CAD development through a number of pathways. It may promote the VSMCs proliferation and migration, promote the injured vascular remodeling, and/or promote the plaque rupture and

However, this meta-analysis is not without limitations. The large-scale studies present in the analysis were not adequate to fully elucidate the complex relationship between the MMP-9-1562C>T gene polymorphism and CAD. The CAD susceptibility is also influenced by environmental factors, such as smoking, diabetes, dyslipidaemia, air pollution, inflammation, and psychological factors (Wang et al., 2012; Agüero et al., 2013; Parruti et al., 2013; Wichmann et al., 2013). It is quite possible that the MMP-9-1562C>T gene polymorphism interact with a risk factor that was not within the scope of this study. There are also many other MMP-9 gene polymorphisms as P574R, R+279Q, and R668Q that influence the MMP-9 serum level (Zhi et al., 2010).

In conclusion, the current meta-analysis indicates that the MMP-9-1562T allele may increase the CAD risk among the Chinese population. This result has the potential to guide the therapy strategy for a CAD patient. Taking into account the This work was funded by the National Natural Science Foundation of China (NSFC 81100073 to YL), Excellent Young and Middle-Aged Teachers Assistance Program of Nanjing Medical University for YL, Jiangsu Overseas Research and Training Program for University Prominent Young and Middleaged Teachers and Presidents, and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). Thank all our colleagues working in the Department of geriatrics, the First Affiliated Hospital of Nanjing Medical University.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2016.00212

## REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Li, Yang, Zhou, Gong, Geng, Kim, Zhou, Qian, Wang and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Stress and Vascular Catastrophes in Newborn Rats: Mechanisms Preceding and Accompanying the Brain Hemorrhages

Oxana Semyachkina-Glushkovskaya1, 2 \*, Ekaterina Borisova<sup>3</sup> , Maxim Abakumov <sup>4</sup> , Dmitry Gorin<sup>5</sup> , Latchezar Avramov <sup>3</sup> , Ivan Fedosov <sup>6</sup> , Anton Namykin<sup>6</sup> , Arkady Abdurashitov <sup>6</sup> , Alexander Serov <sup>1</sup> , Alexey Pavlov <sup>7</sup> , Ekaterina Zinchenko<sup>1</sup> , Vlad Lychagov <sup>6</sup> , Nikita Navolokin<sup>8</sup> , Alexander Shirokov <sup>9</sup> , Galina Maslyakova<sup>8</sup> , Dan Zhu<sup>10</sup> , Qingming Luo<sup>10</sup>, Vladimir Chekhonin<sup>4</sup> , Valery Tuchin2, 6, 11 and Jürgen Kurths 2, 12, 13

#### Edited by:

*Alexey Goltsov, Abertay University, UK*

#### Reviewed by:

*Ningjun Li, Virginia Commonwealth University, USA Denis E. Bragin, University of New Mexico School of Medicine, USA*

#### \*Correspondence:

*Oxana Semyachkina-Glushkovskaya glushkovskaya@mail.ru*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *28 January 2016* Accepted: *22 May 2016* Published: *14 June 2016*

#### Citation:

*Semyachkina-Glushkovskaya O, Borisova E, Abakumov M, Gorin D, Avramov L, Fedosov I, Namykin A, Abdurashitov A, Serov A, Pavlov A, Zinchenko E, Lychagov V, Navolokin N, Shirokov A, Maslyakova G, Zhu D, Luo Q, Chekhonin V, Tuchin V and Kurths J (2016) The Stress and Vascular Catastrophes in Newborn Rats: Mechanisms Preceding and Accompanying the Brain Hemorrhages. Front. Physiol. 7:210. doi: 10.3389/fphys.2016.00210* *<sup>1</sup> Department of Physiology of Human and Animals, Saratov State University, Saratov, Russia, <sup>2</sup> Huazhong University of Science and Technology, Wuhan, China, <sup>3</sup> Laboratory of Biophotonics, Institute of Electronics, Bulgarian Academy of Sciences, Sofia, Bulgaria, <sup>4</sup> Medico-Biological Department, Russian National Research Medical University, Moscow, Russia, <sup>5</sup> Department of Nanotechnology, Saratov State University, Saratov, Russia, <sup>6</sup> Department of Physics, Saratov State University, Saratov, Russia, <sup>7</sup> Department of Electrical Engineering and Electronics, Saratov State Technical University, Saratov, Russia, <sup>8</sup> Department of Pathological Anatomy, Saratov State Medical University, Saratov, Russia, <sup>9</sup> Saratov Research Center, Institute of Biochemistry and Physiology of Plants and Microorganisms, Russian Academy of Sciences (IBPPM RAS), Saratov, Russia, <sup>10</sup> Britton Chance Center for Biomedical Photonics, Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan, China, <sup>11</sup> Laboratory of Biophotonics, Science Department, Tomsk State University, Tomsk, Russia, <sup>12</sup> Department of Physics, Humboldt University, Berlin, Germany, <sup>13</sup> Research Domain Transdisciplinary Concepts and Methods, Potsdam Institute for Climate Impact Research, Potsdam, Germany*

In this study, we analyzed the time-depended scenario of stress response cascade preceding and accompanying brain hemorrhages in newborn rats using an interdisciplinary approach based on: a morphological analysis of brain tissues, coherent-domain optical technologies for visualization of the cerebral blood flow, monitoring of the cerebral oxygenation and the deformability of red blood cells (RBCs). Using a model of stress-induced brain hemorrhages (sound stress, 120 dB, 370 Hz), we studied changes in neonatal brain 2, 4, 6, 8 h after stress (the pre-hemorrhage, latent period) and 24 h after stress (the post-hemorrhage period). We found that latent period of brain hemorrhages is accompanied by gradual pathological changes in systemic, metabolic, and cellular levels of stress. The incidence of brain hemorrhages is characterized by a progression of these changes and the irreversible cell death in the brain areas involved in higher mental functions. These processes are realized via a time-depended reduction of cerebral venous blood flow and oxygenation that was accompanied by an increase in RBCs deformability. The significant depletion of the molecular layer of the prefrontal cortex and the pyramidal neurons, which are crucial for associative learning and attention, is developed as a consequence of homeostasis imbalance. Thus, stress-induced processes preceding and accompanying brain hemorrhages in neonatal period contribute to serious injuries of the brain blood circulation, cerebral metabolic activity and structural elements of cognitive functio These results are an informative platform for further studies of mechanisms underlyin stress-induced brain hemorrhages during the first days of life that will improve the futur generation's health. n. g e

Keywords: stress, newborn rats, cerebrovascular catastrophes, mechanisms

#### INTRODUCTION

Cerebrovascular catastrophes, such as spontaneous brain hemorrhages in apparently normal full term newborns, are a real phenomenon. The problem is that in many cases different types and severity of brain hemorrhages in full term newborns are presented usually without any outward clinical symptoms or with any subtle non-specific neurological signs (Whitby et al., 2004; Looney et al., 2007; Rooks et al., 2008; Gupta et al., 2009). Therefore, the precise incidence and distribution in asymptomatic full term newborns is not known. However, during the last decade, improvement of neuroimaging technique made obvious that asymptomatic brain hemorrhages in full term newborns are frequent (from 26% till 46%) (Looney et al., 2007; Rooks et al., 2008). Notice, the majority of survived full term newborns after brain hemorrhages have a complete recovery (Jhawar et al., 2005; Siu et al., 2006; Rooks et al., 2008; Brouwer et al., 2010). But, neonatal death from brain hemorrhages can reach 25% due to asphyxia (Brouwer et al., 2010). The information about the neurological outcome after brain hemorrhages in full term newborns is extremely limited and typically is focused on the study of relatively short term follow-up with most babies not reaching school-age when higher order deficits are manifest (Jhawar et al., 2005; Siu et al., 2006; Rooks et al., 2008; Brouwer et al., 2010; Takenouchi et al., 2012; Kirton and deVeber, 2013). The results of few studies in this area showed that in future life such babies in 32% of cases develop cognitive deficit (Jhawar et al., 2005) and encephalopathy (Takenouchi et al., 2012), 21%—epilepsy (Siu et al., 2006), 14%—speech delay (Rooks et al., 2008), 8%—cerebral palsy (Brouwer et al., 2010). It is important to note that a child's brain develops over many years. As a result, difficulties can only be detected at an age when it is expected the development of relevant skills. For example, associative learning and thinking might not be recognized for many years. Therefore, the true incidence of cognitive impairment in asymptomatic neonates with brain hemorrhages may be slightly higher than they reported (Semyachkina-Glushkovskaya et al., 2016).

Thus, spontaneous brain bleeding in neonates is a major problem of future generation's health due to the high rate of death and cognitive disability of such newborns. Therefore, the study of mechanisms underlying silent pathological processes preceding and accompanying brain hemorrhages in neonatal period is absolutely essential.

In the majority, the reasons for brain hemorrhages in newborns cannot be found (Gupta et al., 2009). One possible factor is assisted delivery (Benedetti, 1999; Towner et al., 1999) (forceps or vacuum extraction), but those suggestions are not consistent (Whitby et al., 2004; Looney et al., 2007; Rooks et al., 2008). Blood clotting disorders may play an important role in brain bleeding causes in neonatal period (Gupta et al., 2009). However, there are no strong evidences in this detection (Gover et al., 2011). The major risk factor for newborns is stress, which babies have during embryonic development due to different stresses of mother, critical physiological changes during delivery and intensive adaptation to the new conditions of life. The first 3 days of life are the most critical. Mortality is highest in the initial 24 h after birth, up to 50% die within the first 3 days of life, and about 75% of all neonatal deaths occur in the first week of life (early neonatal death; Paul, 2006; Lawn et al., 2010). It is well-known that the neonatal environment is an important determinant of stress-related diseases (Maccari et al., 2003; Mirescu et al., 2004). However, stress itself and mechanisms underlying stress-related brain injury in neonates are not well-studied due to methodological difficulties.

Aiming to achieve a better understanding of mechanisms underlying stress-induced brain hemorrhages in neonatal period, in this experimental study on newborn rats, we intend to uncover a time-depended scenario of stress response cascade preceding and accompanying brain hemorrhages using interdisciplinary approach based on a morphological analysis of brain tissues, coherent-domain optical technologies for visualization of the cerebral blood flow (CBF), monitoring of the cerebral oxygenation and deformability of red blood cells (RBCs).

### METHODS

#### Subjects

Experiments were carried out in newborn mongrel rats, 12 days old using three groups: (1) intact, unstressed newborn rats in the control group (n = 25); (2) stressed rats in the pre-hemorrhage groups (2, 4, 6, 8 h after stress, n = 25 in each group); (3) stressed rats in the post-hemorrhage group (24 h after stress, n = 27). All procedures were performed in accordance with the "Guide for the Care and Use of Laboratory Animals." The experimental protocol was approved by the Committee for the Care and Use of Laboratory Animals at Saratov State University (Protocol H-147, 17.04.2001).

Note, that the rat is a good animal subject for the study of the development of cerebrovascular catastrophes during the first days of life due to similar dynamics of the brain maturation in humans (Coyle, 1977). Our choice of the age of newborn rats is caused by the fact that at postnatal 12 days the brain of rats is close to the development of the brain in full-term human neonates (Rice et al., 1981).

#### Modeling of Brain Hemorrhages

To induce brain hemorrhages, the following protocol of sound stress's impact was used (120 dB, 370 Hz): 10 s of sound followed by a 60-s pause; this cycle was repeated throughout a 2 h period (Semyachkina-Glushkovskaya et al., 2013, 2015a,b).

To analyze the brain-injures induced by stress and to confirm the development of brain hemorrhages, all newborn rats were decapitated for a histological study of brain tissue. The samples were fixed in 10% buffered neutral formalin. The formalin fixed specimens were embedded in paraffin, sectioned (4 µm), and stained with haematoxylin and eosin.

# State-of-the-Art Interdisciplinary Approaches

To study the time-depended scenario of stress response cascade preceding and accompanying brain hemorrhages in newborn rats we used an interdisciplinary approach based on a morphological analysis of brain tissues, coherent-domain optical technologies for visualization of the cerebral blood flow, monitoring of the cerebral oxygenation, and the deformability of red blood cells (RBCs; **Figure 1**).

## Measurement of Cerebral Blood Flow

To assess stress-related changes in cerebral circulation, we used a commercial swept source Doppler optical coherence tomography (DOCT) system OCS1300SS (Thorlabs Inc. USA) operating at 1325 nm central wavelength and 100 nm bandwidth. Transverse and axial resolutions of the DOCT system are 25 and 12 µm (on the air) respectively. A-scan rate is equal to 16 kHz, which allows us to measure absolute velocities up to ∼5.5 mm/s (You et al., 2014).

#### Analysis of Cerebral Oxygen Saturation

The level of blood oxygen saturation (SpO2) in the brain as an important criterion of cerebral metabolic activity of the different functional states of an organism (Liu et al., 2015) was monitored by using the pulse oximeter model CMS60D (Contec Medical Systems Co., Ltd., Qinhuangdao, China). The optical sensor was based on dual wavelengths pulse oximetry approach, using 660 and 880 nm for the SpO<sup>2</sup> detection. The oxy-hemoglobin saturation (SpO2) is given as a percentage of HbO<sup>2</sup> vs. the total Hb in the blood. For confirmation of oximetry results, we used magnetic resonance imaging (MRI, Clin Scan 7T, Bruker Biospin) in susceptibility weighted regime (SWI) for determination of oxygenation of the brain tissues in the pre-hemorrhage (4 h after stress, n = 10) and the post-hemorrhage (24 h after stress, n = 10) groups in comparison with the control group (n = 10). SWI is a novel technique that is highly sensitive to local field in homogeneities and venous deoxygenated blood (Tsui et al., 2009).

# Evaluation of Red Blood Cell Deformability

Micropipette aspiration method was used to measure the resistance to deformation or stiffness of the membranes of red blood cells. It is based on measurement of RBC aspiration depth at given negative pressure applied to a micropipette with internal diameter of 2–2.5 µm (Mitchison and Swann, 1954; Rand and Burton, 1964; Artmann et al., 1997; Sinha et al., 2015; Zheng et al, 2015). Homemade glass micropipettes with 2 µm internal diameter were used for measurements. Pipettes were pulled using 0.93 mm borosilicate glass capillaries washed with isopropanol

and distilled water. Aspiration procedure was imaged with inverted microscope using oil immersion lens (100 × NA=1.25) and color CMOS camera (DCC1615C, Throlabs, Germany) to acquire RBC image. Pressure was applied using U manometer filled buffer solution and connected to the pipette. Blood samples for micropipette aspiration were diluted in 400 times with buffer solution. Aspiration depth 1X was measured in equilibrium conditions. The goal of our study was to assess resistance to deformation for maximal possible number of cells in each sample. To reduce time required for individual measurement we determined only the deformation (delta x) of cell when constant low pressure is applied. After the negative pressure was applied to cell the aspirated part of membrane quickly reaches the equilibrium state when it no longer shifts along the pipette. Delta X was measured after membrane displacement rate decreases below one pixel (50 nm) for 30 s. It enables fast (2–3 min per cell) measurements of 30 cell within 1.5 h after blood sample preparation. That is sufficiently faster than measurements of critical pressure needed for estimation of membrane tension (Sinha et al., 2015).

#### Statistical Analysis

The results were presented as mean ± standard error of the mean (SEM). Differences from the initial level in the same group were evaluated by the Wilcoxon test. Intergroup differences were evaluated using the Mann-Whitney test and ANOVA-2 (post-hoc analysis with the Duncan's rank test). The significance levels were set at p < 0.05 for all analyses.

## RESULTS

# Changes in the Cerebral Venous Hemodynamics Preceding and Accompanying Stress-Induced Brain Hemorrhages in Newborn Rats

The vascular system is one of the first, which responds to stress to provide an adequate metabolism during the mobilization of the organism. However, in the neonatal period the vascular sensitivity to stress is different between perfusion (microcirculatory) and capacitive (venous) sectors of the cerebral circulation. In our previous study on newborn rats, we showed that the stress-reactivity of cerebral veins is higher than that of microvessels due to the specific function of veins to immediately change blood volume during stress and to maintain filling pressure to the heart (Semyachkina-Glushkovskaya et al., 2015a). Here we focused on the study of dynamic changes in the venous component of the cerebral circulation preceding and accompanying stress-induced hemorrhages in newborn rats. For this purpose, we selected the sagittal sinus that is one of major sinuses collecting blood from the small veins of the brain and directs it into the peripheral circulation. We choose this vessel to perform measurements of hemodynamic parameters through the anterior fontanel as a window to the brain in newborn animals.

**Figure 2** shows the DOCT results, which demonstrate timedepended changes in diameter of the sagittal sinus and blood flow velocity before and after stress-induced brain hemorrhages in newborn rats. During the pre-hemorrhage period (2, 4, 6, and 8 h after stress), we observed a increase in the size of the sagittal sinus and a gradual decrease in the blood flow velocity. The maximal changes of these hemodynamic parameters were observed during the incidence of brain hemorrhages (24 h after stress).

#### Changes in the Blood Oxygen Saturation of the Brain Tissues and Erythrocyte Deformability Preceding and Accompanying Stress-Induced Brain Hemorrhages in Newborn Rats

Although investigators have implicated hypoxia as a potential risk factor for the brain hemorrhages in neonates (Michoulas et al., 2011; Luo et al., 2014; van der Aa et al., 2014), the role of hypoxia in the intracranial hemorrhages remains controversial because brain bleeding itself may cause respiratory distress. Therefore, it is difficult to ascertain, if cerebral hypoxia is a key factor in brain hemorrhages, or if it is a consequence of this condition (Jhawar et al., 2003). For a better understanding of the role of hypoxia in the development of brain hemorrhages in the neonatal period, at the second step of our work, we studied the level of SpO<sup>2</sup> in the brain of newborn rats. With this aim, we used the pulse oximetry that confirmed by MRI-SWI imaging.

**Figure 3** shows that the pre-hemorrhage group (n = 15) demonstrated a gradual reduction of the SpO<sup>2</sup> level. Indeed, the decrease in SpO<sup>2</sup> was observed during the whole time before the incidence of brain hemorrhages with maximal reduction of oxygen delivery 4–8 h after stress (by 27, 28, 30, respectively; p < 0.05). The post-hemorrhage group showed the decrease in SpO<sup>2</sup> level by 27% (p < 0.05), i.e., it remained low.

**Figure 4** demonstrates the example of MRI-results in SWI regime obtained from 10 newborn rats in each group (the control, pre- and post-hemorrhage). The pre- and post-hemorrhage periods were characterized by hypointense contrasting of cerebral

FIGURE 3 | Time-dependent changes in blood oxygen saturation of the brain of newborn rats in the pre-hemorrhage (2, 4, 6, 8 h after stress, n = 25 in each group) and post- hemorrhage (24 h after stress, n = 27) groups in comparison with control group (n = 25). \*\*\**P* < 0.001 compared with the control values.

vessels. These changes can be caused by the accumulation of deoxyhemoglobin due to prolonged hypoxia leading to reduced oxygen levels.

The diminished oxygen delivery to the brain tissues in response to a hypoxia may be related to changes of RBC mechanical properties (Ellsworth et al., 1995; Sprague et al., 1996, 2001, 2003; Parthasarathi and Lipowsky, 1999; Kirby et al., 2012). However, the information about the role of erythrocyte deformability in hypoxia after brain hemorrhages is significantly limited (Pollock and Harrison, 1982; Grotta et al., 1986). The question about the pre- and post-hemorrhage changes in the RBCs elasticity during hypoxia remains open. We assumed that the changes of RBC deformability might be one of mechanisms responsible for hypoxia in stressed newborn rats. To check our hypothesis we analyzed RBC deformability in newborn rats before and after stress-induced hemorrhages using the method of micropipette aspiration. We chose two important points for experiment: 2 h after stress when the first decrease in SpO2 occurred in the pre-hemorrhage group (n = 15) and 24 h after stress in the post-hemorrhage group (n = 15).

**Figure 5** demonstrates that the RBCs deformability increased by 10% (p < 0.05) in the pre-hemorrhage group and by 40% (p < 0.001) in the post-hemorrhage group.

FIGURE 5 | Time-dependent changes in the aspiration depth of red blood cells in the pre-hemorrhage (2 h after stress, n = 25) and posthemorrhage (24 h after stress, n = 27) groups in comparison with control group (n = 25). Figure shows aspiration depth at 30 N/m<sup>2</sup> pressure measured in blood samples taken in different time interval after stress. Each point corresponds to a value averaged over measurements performed with 30 different RBCs. \**P* < 0.05; \*\*\**P* < 0.001 compared with the control values.

## Changes in the Brain Areas Involved in Higher Cognitive Functions before and after the Stress-Induced Brain Hemorrhages in Newborn Rats

The prefrontal cerebral cortex is a special zone of the brain, which includes the higher-order association areas and plays a key role in memory, attention, perception, awareness, thought, language, and consciousness. The cognitive deficit in babies is highly associated with brain hemorrhages, which they had in neonatal period in the frontal lobe (Jhawar et al., 2005; Rooks et al., 2008; Semyachkina-Glushkovskaya et al., 2016). However, the mechanisms underlying these pathological processes remain poorly understood and request detailed studies in this field. With this aim, at the final step of our work, we studied the morphological changes in the molecular layer of the prefrontal cortex and in the pyramidal neurons, which are crucial for the "feedback" interactions in the cerebral cortex involved in associative learning and attention (Gilbert and Sigman, 2007).

Our results clearly show that the pre-hemorrhage group (n = 15) demonstrated a gradual decrease in the thickness of the molecular layer of the cortex and in the number of pyramidal neurons with reducing of their diameter (**Figures 6**–**9**). So, during the pre-hemorrhage period (2, 4, 6, and 8 h after stress) the thickness of the molecular layer of the cortex was decreased by 31, 40, 62, and 68%, respectively, p < 0.05; the number of pyramidal cells—by 17, 21, 21, and 33%, respectively, p < 0.05; the diameter of pyramidal neurons—by 19, 20, 29, and 39%, p < 0.05, respectively (**Figures 8**, **9**). The incidence of the brain bleeding was accompanied by a more pronounced decrease in the thickness of the molecular layer of the cortex, which was reduced by 73% (p < 0.05) compared with the normal state. The number of pyramidal cells and their diameter were decreased

in each group) and post- hemorrhage (24 h after stress, n = 27) groups in comparison with control group (n = 25). Hematoxylin and Eosin staining. Bars represent 10 µm (774X). (A)—the normal state (the control); the pre-hemorrhage period (B–E—2, 4, 6, 8 h after stress); the post-hemorrhage period (F—24 h after stress).

by 34% (p < 0.05) and 38% (p < 0.05), respectively, i.e., it remained reduced. These pathological changes in the cortex were accompanied by the apoptosis of pyramidal cells (**Figure 10**). Furthermore, we observed the development of the perivascular edema during all pre- and post-hemorrhages time in all newborn rats (**Figure 11**). But, there were no any stressful changes in severity of the fluid pathway from the cerebral vessels in the preand post-hemorrhage groups.

The results of three series of experiments suggest that the reduction of blood circulation in the cerebral venous system and hypoxia are the important factors for disorders of the functional platform for the integration of brain centers, such as the molecular layer of the cortex and the pyramidal neurons.

## DISCUSSION

In this interdisciplinary study we analyzed the time-dependent scenario of stress-induced brain hemorrhages in newborn rats on the different levels of cascade of stress reaction on morphological (histological analysis of cerebral cortex), systemic (monitoring of cerebral venous blood flow), metabolic (assessment of oxygen saturation of the brain tissues), and cellular (changes in RBCs deformability and pyramidal neuron number and diameter) levels.

In our previous studies, we clearly showed the location, types, depth, and morphological changes related to stress-induced brain hemorrhages in newborn rats (Semyachkina-Glushkovskaya et al., 2015a,b). Based on our earlier studies, we suggest that the venous component of the cerebral circulation is highly sensitive to the stress in the neonatal period (Semyachkina-Glushkovskaya et al., 2013, 2015a; Pavlov et al., 2014). Here we focused on the study of time dynamics of changes in the cerebral venous system preceding and accompanying stressinduced brain hemorrhages in newborn rats. The object of our study was the sagittal sinus, which collects blood from all veins of the brain and directs it into the peripheral circulation. DOCT imaging uncovers time-dependent progressive changes in this vessel during the pre-hemorrhage time (2, 4, 6, and 8 h after stress) and after incidence of brain bleeding (24 h after stress). The latent period of brain hemorrhages is characterized by a relaxation of the sagittal sinus with a fall of blood flow velocity in it that was more pronounced after incidence of brain bleedings in newborn rats. Notice, the changes in these two hemodynamic parameters are closely not correlated. We observed the gradual reduction of venous blood flow during all pre- and post-hemorrhage time, while the increase in the diameter of sagittal sinus was from 6 her till 8 her after stress stable. But, 24 h after stress the sagittal sinus dilated significantly compared with other periods of observations. These changes can be explained by mechanisms underlying distributions of cerebral blood flow under stress. In our previous work we showed that the progressive relaxation of cerebral veins in newborn rats with the stroke causes accumulation of blood not only in venous network but also in microvessels due to the redistribution of blood flow in the cerebral vessels to decrease the pressure of accumulated blood on the thick walls of cerebral veins (Semyachkina-Glushkovskaya et al., 2015a).

We also observed that the pathological relaxation of the main cerebral vein was accompanied by the formation of perivascular edema i.e., fluid pathway from the vessels. In our previous work, the histological data show that the increase in size of the sagittal sinus is associated with the dilation of all cerebral veins, especially in the pail matter of the cerebral cortex (Semyachkina-Glushkovskaya et al., 2015a). The relaxation of cerebral veins with

perivascular edema is a marker of accumulation of an extensive blood in the venous system and suppression of blood outflow from the brain leading to venous insufficiency (Valdueza et al., 2013).

Thus, the sagittal sinus shows sensitive changes to the deleterious effect of stress from the pre-hemorrhage time until the incidence of brain bleeding. Clinical studies also have shown that neonatal intracranial hemorrhages are primary venous infarction due to a weakness of the wall of cerebral veins in neonates (Hambleton and Wigglesworth, 1976; Ghazi-Birry et al., 1997; Bruno et al., 2014).

A gradual reducing of the blood flow velocity in the dilated sagittal sinus in stressed newborn rats is associated with a time-dependent reduction of oxygen supply to the brain. The results of our experiments clearly show that cerebral hypoxia precedes stress-induced brain hemorrhages. This fact is consistent with other experimental and clinical data. Thoresen et al. in experiments on newborn pigs demonstrated that severe hypoxia itself can induce spontaneous brain hemorrhage in newborn pigs (Thoresen et al., 2001). Aderliesten et al. in clinical observations showed a lower cerebral fraction tissue oxygen extraction in newborns before severe stroke (Alderliesten et al., 2013). Taking into account these facts we suppose that hypoxia in the pre-hemorrhage period might be a causative factor provoking critical changes in the brain, associated with intracerebral hemorrhage in newborn rats.

neurons in newborn rats in the pre-hemorrhage (2, 4, 6, 8 h after stress, n = 25 in each group) and post- hemorrhage (24 h after stress, n = 27) groups in comparison with control group (n = 25). \**P* < 0.05; \*\**P* < 0.01; \*\*\**P* < 0.001 compared with the control values.

Our results clearly exhibit that the cerebral hypoxia and disturbances in the cerebral venous system in newborn rats are accompanied by an increase in deformability of RBCs. The deformation of RBCc is a key factor for a release of powerful vasorelaxant such as adenosine 5′ triphosphate (ATP) (Sprague et al., 1996). This mechanism is related to the activation by β-adrenergic receptors presented at the RBCs' surface (Olearczyk et al., 2001). In our previous study, we showed the high sensitivity of cerebral vessels of newborn rats to a pharmacological modulation of vascular β-adrenergic receptors (Pavlov et al., 2014). These facts allow us to believe that the increased deformability of RBCs in stressed newborn

FIGURE 10 | The example of apoptotic bodies (arrowed) in the prefrontal cortex obtained from 27 newborn rat 24 h after stress (the signs of apoptosis are condensation and fragmentation of nuclei in an intensely eosinophilic cytoplasm). Hematoxylin and Eosin staining. Bars represent 10 µm (774X).

from the cerebral vessels (arrowed) in the brain obtained from the pre-hemorrhage group (2, 4, 6, 8 h after stress, n = 25 in each group) and in the post-hemorrhage group (24 h after stress, n = 27).

rats might be one of possible mechanisms contributing to the pathological relaxation of the cerebral veins due to ATP release.

A stress-induced gradual reduction of oxygen supply to the brain and blood flow velocity in the cerebral venous system were accompanied by gradual pathological changes in the brain areas involved in higher cognitive functions such as the molecular layer of the prefrontal cortex and the pyramidal neurons, which are crucial for associative learning and attention (Gilbert and Sigman, 2007). The severity of stressrelated disorders on the level of cerebral venous circulation and oxygenation of the brain tissues was associated with the time-dependent decrease in the thickness of the molecular layer of the prefrontal cortex as well as with the reduction of the number of pyramidal neurons and their diameter. The formation of perivascular edema, which we observed in the pre- and post-hemorrhage periods, can be one of the reasons responsible for a decrease in the diameter of the pyramidal neurons due to the fluid pathway from the cerebral vessels and mechanical compression. Notice that the incidence of the brain hemorrhages was accompanied by the progression of pathomorphological changes in the cerebral cortex up to the apoptosis of pyramidal neurons, i.e., irreversible changes leading to the death of cerebral cells. In our previous morphological studies of the brain tissues and cerebral vessels in newborn rats we clearly show the progressive accumulation of blood in the superficial cerebral veins (the pre- hemorrhage time) and in deep cerebral veins and microvessels (the post-hemorrhage time) (Semyachkina-Glushkovskaya et al., 2015a). We assume that vascular component of stress-related brain injures play an important role in drastic shrinkage of the prefrontal cortex via mechanisms of compressive pressure on the molecular layer of the prefrontal cortex.

What does the brain hemorrhages mean for intellectual functions of the brain? The results of our study give clear evidence that even an early latent stage of brain bleeding is associated with significant morphological lesions of the "intellectual zone" of the prefrontal cortex that is accompanied by the irreversible apoptosis process after the incidence of brain hemorrhages. Thus, stress-induced processes preceding and accompanying brain hemorrhages in neonatal period contribute the serious injures of the brain association areas responsible for cognitive functions. Future studies need to focus on long-term outcomes after brain hemorrhages induced by stress in neonatal period to develop effective prognostic criteria and to optimize neuroprotective strategies.

# CONCLUSION

In general, our results suggest significant time-dependent stressinduced changes in the brain preceding and accompanying intracranial hemorrhages in newborn rats. The pre-hemorrhage processes are realized at different levels of stress response cascade such as systemic (the progressive reduction of blood flow circulation in the cerebral venous system), metabolic (the cerebral oxygenation declines), and cellular one's (the increase in RBCs deformability; the decrease in the thickness of the molecular layer in the prefrontal cortex and in the number of pyramidal neurons with reducing their diameter). The posthemorrhage time is characterized by progression of stressinduced systemic, metabolic, and cellular changes in the brain, which are accompanied by irreversible cell death apoptosis process in the brain areas involved in higher cognitive functions. Thus, the stress-induced processes preceding and accompanying brain hemorrhages in the neonatal period contribute to serious injures of the brain blood circulation, cerebral metabolic activity and the structural elements of cognitive function. These results are an informative platform for further studies of mechanisms underlying stress-induced brain hemorrhages during the first days of life that will improve the health of the future generation.

# AUTHOR CONTRIBUTIONS

OS: wrote the text of article, analyzed all results. EB: performed a measurement of cerebral oxygenation. MA: performed monitoring of oxygen saturation using magnetic resonance imaging (MRI) of the brain of newborn rats. DG: prepared protocol of anesthesia for newborn rats during all experiments. LA: made analysis of the changes in the cerebral oxygenation. IF: developed the plan of optical experiments with the red blood cells deformability. AN: performed a measurement of red blood cells deformability. AA: performed a measurement of cerebral venous blood flow. AS: performed a model of stress-induced brain bleeding in newborn rats. AP: performed the statistical analysis of data. EZ: prepared animals for all experimental tasks. VL: analyzed the recording of DOCT data. NN: performed the histological analysis. GM: analyzed the histological data. DZ: gave consultation about a measurement of cerebral blood flow in newborn animals. QL: discussed DOCT results. VC analyzed MR data. VT: coordinated experiments with the coherent-domain optical technologies. JK: discussed results and application of interdisciplinary methods in the experiments. AS made algorithm of sound stress and selected suitable frequence of sound to induce brain hemorrhages in newborn rats. He also performed the final approval of article.

### ACKNOWLEDGMENTS

The work of OS, IF, AN, AA, AS, AP, EZ, NN, VT, and JK in preparation of this original article, discussion of obtained results, all experimental works including magnetic resonance imaging, measuring of cerebral blood flow using Doppler coherent tomography, oximetry, histological analysis, investigations of RBCs deformability was supported by Grant of Russian Science Foundation no. 14-15-00128. The work of EB, LA in the study of performing of pulse oximetry was supported by Grant DFNI-B02/9/2014 of Bulgarian National Science Fund. The work of DZ, QL related to discussion was supported by National Nature Science Foundation of China under Grants 81171376 and 91232710, the Science Fund for Creative Research Group of China under Grant 61421064, and the Program of Introducing Talents of Discipline to Universities in China under Grant B07038. MRI imaging was carried out on ClinScan 7T machine included into Center for Shared Usage "Medical and biotechnological nanotechnologies."

# REFERENCES


glucocorticoid hormones. Neurosci. Biobehav. 27, 119–127. doi: 10.1016/S0149- 7634(03)00014-9


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Semyachkina-Glushkovskaya, Borisova, Abakumov, Gorin, Avramov, Fedosov, Namykin, Abdurashitov, Serov, Pavlov, Zinchenko, Lychagov, Navolokin, Shirokov, Maslyakova, Zhu, Luo, Chekhonin, Tuchin and Kurths. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Decomposition of Near-Infrared Spectroscopy Signals Using Oblique Subspace Projections: Applications in Brain Hemodynamic Monitoring

Alexander Caicedo1, 2 \*, Carolina Varon1, 2, Borbala Hunyadi 1, 2, Maria Papademetriou<sup>3</sup> , Ilias Tachtsidis <sup>3</sup> and Sabine Van Huffel 1, 2

*<sup>1</sup> Department of Electrical Engineering ESAT, STADIUS Center for Dynamical Systems, Signal Processing, and Data Analytics, KU Leuven, Leuven, Belgium, <sup>2</sup> (iMinds Medical) Department Medical Information Technologies, Leuven, Belgium, <sup>3</sup> Biomedical Optics Research Laboratory, Department of Medical Physics and Bioengineering, University College London, London, England*

#### Edited by:

*Alexey Goltsov, Abertay University, UK*

#### Reviewed by:

*Keum-Shik Hong, University of Illinois at Urbana-Champaign, USA Victor Chernomordik, National Institutes of Health, USA Oxana Semyachkina-Glushkovskaya, Saratov State University, Russia Igor Meglinski, University of Oulu, Finland*

> \*Correspondence: *Alexander Caicedo acaicedo@esat.kuleuven.be*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *25 May 2016* Accepted: *19 October 2016* Published: *08 November 2016*

#### Citation:

*Caicedo A, Varon C, Hunyadi B, Papademetriou M, Tachtsidis I and Van Huffel S (2016) Decomposition of Near-Infrared Spectroscopy Signals Using Oblique Subspace Projections: Applications in Brain Hemodynamic Monitoring. Front. Physiol. 7:515. doi: 10.3389/fphys.2016.00515* Clinical data is comprised by a large number of synchronously collected biomedical signals that are measured at different locations. Deciphering the interrelationships of these signals can yield important information about their dependence providing some useful clinical diagnostic data. For instance, by computing the coupling between Near-Infrared Spectroscopy signals (NIRS) and systemic variables the status of the hemodynamic regulation mechanisms can be assessed. In this paper we introduce an algorithm for the decomposition of NIRS signals into additive components. The algorithm, SIgnal DEcomposition base on Obliques Subspace Projections (SIDE-ObSP), assumes that the measured NIRS signal is a linear combination of the systemic measurements, following the linear regression model y = Ax + ǫ. SIDE-ObSP decomposes the output such that, each component in the decomposition represents the sole linear influence of one corresponding regressor variable. This decomposition scheme aims at providing a better understanding of the relation between NIRS and systemic variables, and to provide a framework for the clinical interpretation of regression algorithms, thereby, facilitating their introduction into clinical practice. SIDE-ObSP combines oblique subspace projections (ObSP) with the structure of a mean average system in order to define adequate signal subspaces. To guarantee smoothness in the estimated regression parameters, as observed in normal physiological processes, we impose a Tikhonov regularization using a matrix differential operator. We evaluate the performance of SIDE-ObSP by using a synthetic dataset, and present two case studies in the field of cerebral hemodynamics monitoring using NIRS. In addition, we compare the performance of this method with other system identification techniques. In the first case study data from 20 neonates during the first 3 days of life was used, here SIDE-ObSP decoupled the influence of changes in arterial oxygen saturation from the NIRS measurements, facilitating the use of NIRS as a surrogate measure for cerebral blood flow (CBF). The second case study used data from a 3-years old infant under Extra Corporeal Membrane Oxygenation (ECMO), here SIDE-ObSP decomposed cerebral/peripheral tissue oxygenation, as a sum of the partial contributions from different systemic variables, facilitating the comparison between the effects of each systemic variable on the cerebral/peripheral hemodynamics.

Keywords: oblique subspace projections, Tikhonov regularization, biomedical signal processing, NIRS, brain hemodynamics

## 1. INTRODUCTION

Signal decomposition methods aim at representing a signal as a combination of components that fulfill some specific criteria. For instance, a wavelet transform decomposes a signal into a set of time series that are localized both in frequency and time. Generally, these decomposition schemes define a subspace using an orthonormal basis onto which the signal to be decomposed is projected. In this way the signal is represented in the transformed subspace as a linear combination of the given basis vectors. These decomposition schemes can be seen as a linear regression problem of the form **y** = **Ax** + ǫ, where **y** represents the signal of interest, **A** is a matrix containing the basis for the subspace defined by the transformation, **x** is a vector containing the decomposition coefficients that represent the signal in the transformed domain, and ǫ represents the error term.

Signal decomposition methods can also be linked to the field of system identification, where the output of a system is computed as a linear combination of the partial contributions of the input variables (Ljung, 1999). An identification problem can be formulated as a regression problem of the same form mentioned above, **y** = **Ax** + ǫ, which is a matrix representation of a convolution operation, with **x** being the impulse response of the system. However, in this case, the matrix **A** is constructed differently, and does not contain an orthonormal basis representing the transformed subspace. Instead, this matrix is created from input-output observations of the system using a specific model form to fit the output variable. Models such as ARX (autoregressive with exogenous input), ARMAX (autoregressive moving average with exogenous input), ARIMA (autoregressive integrative moving average), among others have been extensively studied in the literature (Ljung, 1999). These models can be used for simulation, as well as prediction of the system output **y**, e.g., in the prediction of the brain hemodynamic response to an stimulus (Aqil et al., 2012). Others examples of the use of this models in the field of near infrared spectroscopy can be found in Hong and Naseer (2016), Kamran and Hong (2013), Pillonetto (2016), Kamran and Hong (2014), Naseer and Hong (2015). However, the main focus of these models is not to compute the individual contribution of each input variable on the output, but to produce a good estimation of the behavior of the system as a whole, which leads to a lack of interpretability of the models in terms of the underlying physiological processes.

Since real-life problems are likely to be characterized by correlated inputs, the system subspace will consist of a set of nonorthogonal subspaces, challenging the identification of individual contributions, as well as obscuring its clinical interpretation. For instance, when using least squares to identify the contribution of a single subsystem on the output, the noise is assumed to be orthogonal to the system subspace, which due to the presence of correlated inputs is clearly not the case, hence producing erroneous estimates. This problem can be mitigated by the use of an appropriate projector that allows a more effective separation of the different subsystems' dynamics, such as an oblique subspace projector.

Oblique subspace projectors (ObSP) are projection matrices that use a reference subspace trajectory to project a signal onto a desired subspace. When the reference subspace is orthogonal to the desired subspace ObSP reduces to an orthogonal projector. ObSP possesses properties that can be exploited to produce an appropriate decomposition of biomedical signals, in terms of the linear contributions of the different input variables. This allows to provide physical interpretation to the system, in the framework of identification models.

Specifically, this kind of decomposition framework can be used to decipher the relation between different biomedical signals, helping to understand their relations, as well as facilitating their clinical interpretation. In particular in the current clinical practice where patients measurements comprise a large number of concomitant measurements of different biomedical signals, ObSP can aid in facilitating a more accurate interpretation of physiological or patho-physiological processes, since the influence of each input variable on the output can be analyzed separately without interference from confounding factors, this is a critical factor for the clinical interpretability of mathematical models (Slinker and Glantz, 1985). For instance, for the assessment of cerebral autoregulation in premature infants surrogate measurements of cerebral blood flow (CBF), obtained using Near-infrared Spectroscopy (NIRS), can be used. But only when there are not strong variations in SaO<sup>2</sup> (Soul et al., 2007; Wong et al., 2008). However, NIRS measurements are highly coupled to changes in arterial oxygen saturation (SaO2), which introduces information that is not directly linked to changes in CBF. This hinders the use of NIRS as a technology for the bed-side assessment of cerebral autoregulation. In order to correctly assess the status of the cerebral autoregulation mechanism, using NIRS, changes in SaO<sup>2</sup> should be decoupled from the NIRS measurements, prior to further processing. In such an example the use of a decomposition framework, as the one presented in this paper, becomes relevant since it can be used as a preprocessing step to decouple the influences of changes in SaO<sup>2</sup> from the NIRS signals.

Applications of ObSP are scarce in signal processing. Among the applications we can find in the literature we list the work from Behrens and Scharf (1994), who presents the use of ObSP for interpolation, decoding, and elimination of symbol interferences in a communication channel. Tu et al. (1999), used ObSP for hyperspectral image classification. They applied ObSP to quantify the mixture of spectral signatures, from different materials, contained in a specific pixel of an hyperspectral image. Van Overschee and De Moor (1994) proposed the used of subspace system identification, which intrinsically uses ObSP in order to estimate the state space of the system. In the biomedical signal processing field, in some previous work, we have proposed the use of ObSP in combination with wavelet decomposition for cerebral hemodynamics monitoring (Caicedo et al., 2012). There, we aimed at decomposing cerebral hemodynamic signals, measured by means of NIRS, as a sum of the partial linear contributions of different systemic variables such as, mean arterial blood pressure (MABP), SaO2, heart rate (HR), end tidal CO2, among others. We also proposed the use of ObSP as a preprocessing method for NIRS measurements (Caicedo et al., 2013a), and for the extraction of features in a sleep apnea detection algorithm (Varon et al., 2015). In this paper we introduce the theoretical framework for a decomposition algorithm based on ObSP which can be applied in the field of cerebral hemodynamic monitoring. In addition we impose smoothness in the regression parameters by the use of Tikhonov regularization using a matrix differential operator, in contrast with the wavelet decomposition framework of our previous work, resulting in a more formal and flexible problem definition, a more robust solution, as well as a more clear decomposition framework. We test the performance of ObSP with a synthetic example as well as with 2 application examples in the biomedical field, specifically related to the monitoring of brain hemodynamics regulation. In the first application example we decouple the changes in SaO<sup>2</sup> from the NIRS recordings, the SaO<sup>2</sup> is measured using a pulse oxymeter. In the second application example, we find the partial linear contributions of each input variable on the changes of cerebral and peripheral tissue oxygenation. In this example we use as input variables concomitant measurements of MABP, EtCO2, HR, SaO2, and ECMO flow.

The rest of the paper is organized as follows. In Section 2 we briefly introduce ObSP. In Section 2.1 we present the geometrical interpretation for ObSP, and we propose to solve the ObSP problem by solving an alternative least squares problem. The new problem definition allows to smoothly introduce Tikhonov regularization. In Section 3 we present the proposed general algorithm for signal decomposition using ObSP. In Section 4 we show the results from the synthetic and applications examples. Finally, in Sections 5, 6 we discuss our main findings and present the concluding remarks.

Along the manuscript we will represent scalars by lowercase variables such as x. Vectors will be represented by bold lowercase variables such as **x**. Matrices will be represented by bold uppercase variables such as **A**. Generic vector subspaces will be represented by blackboard bold letters such as R, and vector subspaces generated by a given matrix will be represented using calligraphy type letter such as V.

#### 2. OBLIQUE SUBSPACE PROJECTIONS

An oblique subspace projection matrix (ObSP) is a linear operator that projects a given vector onto a target subspace following the direction of a reference subspace. Conversely, an orthogonal subspace projection (OrSP) can be seen as a special case of an ObSP where the target and reference subspaces are orthogonal (Yanai et al., 2011), this permits the construction of OrSP using only a basis for the target subspace. In contrast with ObSP, OrSP is the most popular projection operator since it arises naturally in the least squares solution of a linear regression problem.

In general, in order to construct an oblique projector the basis for the target and the reference subspaces should be known. Let V ⊂ R <sup>N</sup> represent the subspace spanned by a matrix **A** ∈ R N×p , where N and p represents the number of rows and columns of the matrix **A** respectively, V<sup>k</sup> ⊂ V the signal subspace spanned by a partition of **A**, **A**<sup>k</sup> , with **A** = [**A**<sup>k</sup> **A**(k) ]. If V = V1⊕V2⊕...⊕Vd, with ⊕ being the direct sum operator, then the oblique projector onto V<sup>k</sup> along V(k) = V<sup>1</sup> ⊕ ... ⊕ V<sup>k</sup> <sup>−</sup> <sup>1</sup> ⊕ V<sup>k</sup> <sup>+</sup> <sup>1</sup> ... ⊕ Vd, denoted by Pk.(k) , with d ≤ p, is given by:

$$\mathbf{P}\_{k,(k)} = \mathbf{A}\_k (\mathbf{A}\_k^T \mathbf{Q}\_{(k)} \mathbf{A}\_k)^\dagger \mathbf{A}\_k^T \mathbf{Q}\_{(k)},\tag{1}$$

where † represents the generalized inverse, d represents the number of signal subspaces embedded in **A** and satisfies d ≤ p, and **Q**(k) is the orthogonal projector onto Null(**A** T (k) ) ⊂ V ⊥ (k) , which is computed as:

$$\mathbf{Q}\_{(k)} = \mathbf{I}\_N - \mathbf{P}\_{(k)},\tag{2}$$

where **P**(k) = **A**(k) (**A** T (k) **A**(k) ) †**A** T (k) is the orthogonal projector onto V(k) ≡ Span(**A**(k) ) (Yanai et al., 2011).

ObSP arises naturally from the solution of a generalized least squares (GLS) estimation problem. The GLS regression problem is defined as:

$$\begin{array}{ll}\min & \mathbf{\epsilon}^T \Sigma \mathbf{\epsilon} \\ \text{s.t.} & \mathbf{\epsilon} = \mathbf{y} - \mathbf{A} \mathbf{x}, \end{array} \tag{3}$$

where 6 = Var[ǫ|**A**], and **A** = [**a** T 1 ;. . .; **a** T N ], where **A** ∈ R N×p , with N equal to the number of observations. Since the cost function is quadratic in terms of the model parameters **x**, the solution has the following closed form:

$$
\hat{\mathbf{x}}\_{GLS} = (\mathbf{A}^T \mathbf{Z}^{-1} \mathbf{A})^\dagger \mathbf{A}^T \mathbf{Z}^{-1} \mathbf{y}.\tag{4}
$$

Substituting Equation (4) in (3) we obtain **y**ˆ = **A**(**A** <sup>T</sup>6 <sup>−</sup>1**A**) †**A** <sup>T</sup>6 −1 **y**, which can be written as **y**ˆ = **Zy**. The matrix **Z** is idempotent, **Z** <sup>2</sup> = **Z**, and not symmetric, **Z** <sup>T</sup> 6= **Z**, which indicates that **Z** is an ObSP matrix, provided 6 <sup>−</sup><sup>1</sup> = **Q**(k) (Yanai et al., 2011).

GLS considers that the noise variance is not constant and mitigates its effects by including a weighting matrix that is obtained from the variance of the estimated noise (Kariya and Kurata, 2004).

The noise has a large influence in the performance of ObSP. Behrens and Scharf (1994), studied the influence of structured and white noise in the projections when using ObSP. They found that, when the signal subspace for the structured noise is known, an ObSP operator can be created such that its kernel contains this subspace, eliminating the influence of the structured noise from the estimation. However, even though ObSP can effectively remove the structured noise, some components from the white noise might be amplified. They concluded that ObSP is more efficient when the structured noise is dominant to the background white noise. Therefore, by reducing the background noise the efficiency of ObSP to remove structured noise can be improved. Normally this would not be a problem when the noise is considered to be another physiological measurement. However, noise can be reflected in small eigenvalues which might inflate the estimation producing results with a high variance. This problem can be mitigated by the use of a regularization term that numerically stabilizes the solution. This approach will be discussed in Section 2.2. First we will discuss the geometrical interpretation for ObSP.

#### 2.1. Geometrical Interpretation

ObSP operators also appear naturally in the least squares solution of a modified linear regression problem. Consider the linear regression problem **y** = **Ax**+ǫ. By multiplying it by **Q**(k) , defined in Equation (2), we obtain the following regression problem:

$$\mathbf{Q}\_{(k)}\mathbf{y} = \mathbf{Q}\_{(k)}\mathbf{A}\mathbf{x} + \mathbf{Q}\_{(k)}\boldsymbol{\epsilon},\tag{5}$$

and its Least Squares solution is then given by:

$$
\hat{\mathbf{x}}\_{LS} = (\mathbf{A}^T \mathbf{Q}\_{(k)} \mathbf{A})^{-1} \mathbf{A}^T \mathbf{Q}\_{(k)} \mathbf{y}. \tag{6}
$$

Replacing Equation (6) in (5), and using the ObSP operator in Equation (1) we obtain:

$$\mathbf{Q}\_{(k)}\mathbf{y} = \mathbf{Q}\_{(k)}\mathbf{P}\_{k,(k)}\mathbf{y} + \mathbf{Q}\_{(k)}\boldsymbol{\epsilon}.\tag{7}$$

Based on Equation (7), we can interpret an ObSP operator as a linear mapping of a given vector, **y**, onto a desired subspace, V<sup>k</sup> ≡ Span{**A**<sup>k</sup> }, following a reference subspace, V(k) ≡ Span{**A**(k) }, such that the mapped vector, **P**k.(k)**y**, has the same orthogonal projection onto the complement of the reference subspace, V ⊥ (k) ≡ Null{**A** T (k) }, as the vector **y**; notice that **A**<sup>k</sup> 6⊂ **A**(k) . This can be clearly seen in **Figure 1**.

Interestingly, the least squares solution to the problem in Equation (5), **x**LS, is the same as the solution to

$$\mathbf{P}\_{k.(k)}\mathbf{y} = \mathbf{P}\_{k.(k)}\mathbf{A}\mathbf{x} + \mathbf{P}\_{k.(k)}\boldsymbol{\epsilon}.\tag{8}$$

This is proven as follows:

PROOF. The Least Squares solution to Equation (5) is equal to **x**ˆLS = ((**Q**(k)**A**) <sup>T</sup>**Q**(k)**A**) −1 (**Q**(k)**A**) T y, since **Q**<sup>T</sup> (k)**Q**(k) = **Q**(k) , **Q**(k)**A** = **Q**(k)**A**<sup>k</sup> , and **Q**<sup>T</sup> (k) = **Q**(k) then the solution simplifies

to **x**ˆLS = (**A** T <sup>k</sup> **Q**(k)**A**<sup>k</sup> ) <sup>−</sup>1**A** T <sup>k</sup> **Q**(k)y. Similarly, solving Equation (8), we obtain **x**ˆLS = ((**P**k.(k)**A**) <sup>T</sup>**P**k.(k)**A**) −1 (**P**k.(k)**A**) <sup>T</sup>**P**k.(k)y. Since **P**k.(k)**A** = **A**<sup>k</sup> , and Pk.(k) is defined in Equation (1), simplifying the equations we obtain **x**ˆLS = (**A** T <sup>k</sup> **Q**(k)**A**<sup>k</sup> ) <sup>−</sup>1**A** T <sup>k</sup> **Q**(k)y.

#### 2.2. ObSP and Tikhonov Regularization

Tikhonov regularization is the most general form of the Ridge regularization, where a linear regression problem is formulated as follows (Golub et al., 1999):

$$\begin{array}{ll}\min\_{\mathbf{x}} & \boldsymbol{\epsilon}^T \boldsymbol{\epsilon} + \mathbf{y} \, \|\Gamma \mathbf{x}\|\_2^2\\ \text{s.t.} & \boldsymbol{\epsilon} = \mathbf{y} - \mathbf{A} \mathbf{x},\end{array} \tag{9}$$

when the matrix Ŵ = **I**, the problem is known as Ridge regression. Tikhonov regularization is used in order to impose some kind of property to the solution vector **x**. In this manuscript the solution vector **x**, due to the construction of the matrix **A**, will represent the impulse response of a physiological system. Therefore, we are interested in imposing smoothness to the estimated output of the linear problem, since it is expected for this response to be smooth. In this context, the regularization matrix Ŵ will take the form of a differential operator, with entries such as Ŵ(i, i) = −1 and Ŵ(i, i + 1) = 1. In order to obtain a regularized output for the oblique projections we can solve the equivalent problem presented in Equation (5):

$$\begin{array}{ll}\min\_{\mathbf{x}} & \hat{\boldsymbol{\epsilon}}^T \hat{\boldsymbol{\epsilon}} + \mathbf{y} \, \|\Gamma \mathbf{x}\|\_2^2\\ \text{s.t.} & \hat{\boldsymbol{\epsilon}} = \mathbf{Q}\_{(k)} \mathbf{y} - \mathbf{Q}\_{(k)} \mathbf{A} \mathbf{x},\end{array} \tag{10}$$

with ǫˆ = **Q**(k)ǫ. The solution, xˆ = ˆxObSP, is given by:

$$\hat{\mathbf{x}}\_{Obs} = (\mathbf{A}\_k^T \mathbf{Q}\_{(k)} \mathbf{A}\_k + \boldsymbol{\nu} \, \mathbf{T}^T \mathbf{F})^{-1} \mathbf{A}\_k^T \mathbf{Q}\_{(k)} \mathbf{y},\tag{11}$$

where **A**<sup>k</sup> represents a partition of **A** containing only the columns that contains the kth regressor, **Q**(k) is computed using Equation (2), γ is the regularization constant, and Ŵ represents the regularization matrix. This problem does not require the computation of ObSP projectors, which are more costly than the construction of orthogonal projectors. In addition, since the norm of oblique projectors might be larger than one, resulting in the possible amplification of noisy components, solving Equation (11) is numerically more stable.

#### 3. SIGNAL DECOMPOSITION BASED ON OBLIQUE SUBSPACE PROJECTIONS (SIDE-OBSP)

Let's consider N observations from a linear time-invariant (LTI), multiple-input and single-output (MISO) system, with output **y** and d input variables **S** = {**s**i} d i=1 , with **s**<sup>i</sup> ∈ R <sup>N</sup>. The output of this system in discrete time can be represented as follows:

$$\mathbf{y}[n] = \sum\_{i=1}^{d} \sum\_{m=0}^{p\_i} \mathbf{s}\_i[n-m]\mathbf{h}\_i[m] + \epsilon,\tag{12}$$

where **h**<sup>i</sup> represents the impulse response of the subsystem that links the input variable **s**<sup>i</sup> and its corresponding contribution to the output **y**<sup>i</sup> , p<sup>i</sup> is the length of the impulse response, and ǫ is the background noise. Notice that the number of input variables, d, also represents the number of signals subspaces spanned by the system. If the system is stable its impulse response, **h**, decays toward zero. Therefore, **h** can be truncated by an appropriate sample number p, such that hi[n] > δ ∀i, with δ being an appropriately chosen threshold, and n = {1, . . . , p}, the selection of p will be discussed later. Then, the model in Equation (12) can be represented as a linear regression problem of the form:

$$\mathbf{y} = \mathbf{A}\mathbf{h} + \boldsymbol{\epsilon},\tag{13}$$

where **y** ∈ R <sup>N</sup> is the output of the system, **A** ∈ R N×(dp) is the regressor matrix that represents the signal subspaces, and **h** ∈ R dp is a vector containing the concatenated impulse responses **h** <sup>T</sup> = {**h** T i } d i = 1 . Using the input matrix **S** = [**s**1, . . . ,**s**d] and approximating **y** as the output of a moving average (MA) system with finite impulse response, we construct the matrix **A** as a block Hankel expansion of the input matrix **S**, using p delayed samples from each input variable, **s**<sup>i</sup> , in order to define its signal subspace.

Given a vector **s**<sup>k</sup> ∈ R <sup>N</sup>, a "Hankel" matrix constructed from **s**<sup>k</sup> consists of forming a square matrix **A**<sup>k</sup> , such that its entries follow **A**<sup>k</sup> (i, j) = **A**<sup>k</sup> (i − 1, j + 1). The columns of the matrix **A**<sup>k</sup> are delayed versions of the vector **s**<sup>k</sup> . In the proposed framework, the matrix **A**<sup>k</sup> is a truncated version where the number of columns per each input signal is limited to the order of the MA model, p. This matrix is called a block Hankel matrix. The order p can be found using cross-validation, or can be set to a sufficiently large number such that all the impulse responses at the sample number p are smaller in magnitude than a selected threshold. The matrix **A** can be obtained by concatenaiting the matrices **A**<sup>k</sup> obtained by the block Hankel expansion of the input signals **s**<sup>k</sup> for k = {1, . . . , d}, such that **A** = [**A**1, . . . , **A**d].

Since we are interested in finding the linear contribution of each input variable **s**<sup>i</sup> on the output **y**, we can reformulate the model in Equation (13) as follows:

$$\mathbf{y} = \mathbf{A}\_k \mathbf{h}\_i + \mathbf{A}\_{(k)} \mathbf{h}\_{(i)} + \boldsymbol{\epsilon},\tag{14}$$

where **A**<sup>k</sup> is the partition of **A** that is related to the input **s**<sup>i</sup> , and it spans the signal subspace of that input. Notice that due to the Hankel expansion of the matrix **S**, the signal subspace corresponding to the ith input variable, **s**<sup>i</sup> , now corresponds to the kth partition of the matrix **A**, **A**<sup>k</sup> ∈ R N×p . **h**<sup>i</sup> ∈ R p represents the impulse response of the system linked to the input **s**i , **A**(k) ∈ R N×(d − 1)p represents the remaining columns of **A**, and **h**(i) ∈ R (d − 1)p is the vector containing the impulse responses of the remaining subsystems. In order to eliminate the influence of the undesired inputs, we can project **y** onto the Span{**A**<sup>k</sup> }. However, the subspace spanned by the residual components **A**(k) = [**A**1, . . . , **A**<sup>k</sup> <sup>−</sup> <sup>1</sup> , **A**<sup>k</sup> <sup>+</sup> <sup>1</sup> , . . . , **A**d] is likely to be oblique to the Span{**A**<sup>k</sup> }. Notice that **A**(k) is obtained by concatenating the block Hankel expansion of all the inputs except the kth input. To solve this problem we can modify it by multiplying Equation (14) by **Q**(k) , computed as in Equation (2). This results in:

$$\begin{array}{l} \mathbf{Q}\_{(k)} \mathbf{y} = \mathbf{Q}\_{(k)} (\mathbf{A}\_{k} \mathbf{h}\_{k} + \mathbf{A}\_{(k)} \mathbf{h}\_{(k)} + \epsilon), \\ \quad = \mathbf{Q}\_{(k)} \mathbf{A}\_{k} \mathbf{h}\_{k} + \eta\_{k}, \end{array} \tag{15}$$

where **Q**(k)**y** can be seen as a preprocessed version of **y**, where the linear contribution of the regressors **A**(k) has been eliminated, and η<sup>k</sup> = **Q**(k)ǫ is the residual noise component. The solution of this problem using the Tikhonov regularization is given by Equation (11). The linear contribution of the ith input regressor can be found as **<sup>y</sup>**ˆ<sup>i</sup> <sup>=</sup> **<sup>A</sup>**k**h**<sup>ˆ</sup> i . Since this is equivalent to **y**ˆ<sup>i</sup> = **P**k.(k)**y**, and taking into account that oblique projectors might amplify noisy components, as mentioned before, we should guarantee that the structural noise is larger than the background noise to avoid these problems. In the framework of decomposition of biomedical signals, we are interested in removing physiological interferences, hence the structural noise will consist of physiological measurements that normally have a higher power than the background noise.

The decomposition algorithm proposed in this paper, called from now on SIDE-ObSP, is summarized in the Algorithm presented below.

Since the column subspace of the oblique projector Pk.(k) is spanned by the dynamics of the kth input variable, and its Null subspace is represented by the signal subspace from all the other inputs, the oblique projector is able to decouple the dynamics of the kth input from all the other variables, even in the presence of high correlations among them. However, it is important to notice than in the pathological case when the signal subspace of the kth variable is also contained in the Null subspace of the remaining inputs, then the problems becomes ill-defined and the projection matrix tends to infinite. In practice this is counter intuitive, because in this case the Column subspace of the projector will be contained in its Null subspace. In the case of extremely high correlated signals, this might lead to numerical problems that affect the solution, due to the presence of the inverse in Equation

#### **Algorithm. SIgnal DEcomposition based on Oblique Subspace Projections (SIDE-ObSP)**

**Input:** regressor matrix **S** ∈ R N×d , output vector **y** ∈ R N.

**Output:** Matrix formed with the decomposition of the output vector **Y**ˆ = [**y**ˆ1, . . . , **y**ˆd], with **Y**ˆ ∈ R (N − p)×d .


For numerical stability N ≫ dp. A rule of thumb for the number of observations needed is to use 70% of the available data to compute the model parameters and 30% for validation. Taking this into account the number of data points needed for the algorithm, defining an upper limit for p = P, will be N = 10dP, and Nv = 4dP.

(6), however, the use of Tikhonov regularization should be able to partially mitigate this problem as can be seen from Equation (11).

#### 4. APPLICATIONS

#### 4.1. Simulation Study

We considered N = 1024 observations {y<sup>i</sup> , **x**i}, with y<sup>i</sup> ∈ R and **x**<sup>i</sup> ∈ R 3 , following the model **y** = f1(**x**1) + f2(**x**2) + f3(**x**3) + η, where η is uniformly distributed random noise with zero mean, and a chosen standard deviation such that we obtain a signalto-noise ratio SNR = 4db in the signal **y**. The function f<sup>1</sup> was selected as a 3rd order low-pass Butterworth filter, with normalized frequency of 0.15 half-cycles/sample, the functions f2, and f<sup>3</sup> had normalized frequencies of [0.1–0.2] half-cycles/sample and [0.05–0.15] half-cycles/sample, respectively.

We considered the inputs {**x**1, **x**2, **x**3}, to be three binary pseudorandom signals (PBRN) in the normalized frequency range [0–0.3] half-cycles/sample. In addition we induced a correlation of ρ = 0.5 between the input signals by adding a 4th reference PBRN signal, **x**4, to all the inputs using the following formula:

$$\mathbf{x}\_{i} = \rho \mathbf{x}\_{4} + (\sqrt{1 - \rho^{2}})\mathbf{x}\_{i}; \quad i = \{1, 2, 3\}. \tag{16}$$

Before applying SIDE-ObSP, we also contaminated all the inputs with uniformly distributed random noise N(0, σ 2 ), where σ was chosen to reach a SNR = 4db. In addition, to complicate more the decomposition problem, we filtered the signal **x**<sup>4</sup> using a low pass Buttherworth filter of 3rd order and normalized cut-off frequency of 0.3 half-cycles/sample, the filtered signal was mixed together with the output signal **y**, imposing a correlation of 0.3 between them, using Equation (16). We applied ObSP using the noisy inputs, **x**<sup>i</sup> and the noisy output **y** to find the linear contribution of each input on the output, **y**<sup>i</sup> = fi(**x**i).

In **Figure 2**, the estimated impulse responses of the subsystems, **h**ˆ i , using SIDE-ObSP and its unregularized version are shown. We performed 100 simulations and plotted the median and the 95% confidence interval for the estimated impulses responses using regularized and unregularized SIDE-ObSP. It can be seen that regularization smooths the estimations of the impulse response of the systems.

We compared the output from SIDE-ObSP with the output of four different system identification model structures: subspace system identification, ARMAX, ARX and an adaptive filtering model. We used the System Identification toolbox from MATLAB in order to estimate the model. For the subspace system identification we used the function "n4sid" in MATLAB, using zero delays in the input variables and evaluating the order of the system between 2 and 8, since the order of the systems to be identified lie within this range, the final order of the model was selected based on the suggestion given by the algorithm output. For the ARX model, we evaluated different model orders and selected the order that minimizes the Akaike information criteria, as proposed by the function selstruc. For the ARMAX model we fixed the model parameters to [na, n<sup>b</sup> , n<sup>c</sup> , n<sup>k</sup> ] = [5, 5, 5, 0]<sup>1</sup> . We checked that the selected order for the ARMAX models produce a satisfactory output. Finally, for the adaptive filtering approach, we used the function adaptfilt from MATLAB, using the LMS FIR adaptive filter algorithm. The length of the adaptive filter was selected as the identified optimal order of SIDE-ObSP; in addition, the LMS step size was set to 0.008.

To evaluate if SIDE-ObSP is able to effectively decompose the output measurement in a set of linear contributions from each input variable, and if it outperforms the tested system identification models, we computed the estimated output for each one of the three different sub-systems, using the different methods, and computed the cross-entropy between the estimated linear contributions on the output from each input and the different inputs. In short, cross-entropy is a measure of the information transfer between 2 signals. A cross-entropy value equal to 0 indicates no information transfer, cross-entropy values different from zero indicate the amount of information that is transferred from one signal to another, with larger values indicating stronger transfer. For more details on this measure we refer the reader to Faes et al. (2011). We will use from here on

<sup>1</sup>For the interpretation of the variables [na, n<sup>b</sup> , n<sup>c</sup> , n<sup>k</sup> ], please refer to the documentation of the function armax in MATLAB.

the information transfer as a measure for the coupling between the variables. The results are presented in **Table 1**. It can be seen that SIDE-ObSP increases the coupling between each input variable and their corresponding estimated linear influence, whilst reducing the coupling with other input variables. The bold numbers indicate the values of cross-entropy which are different from zeros, i.e., they indicate the pair of signals which present at some degree a linked dynamics. As can be seen in the table only the diagonal elements, in the case of the output provided by SIDE-ObSP, present cross-entropy values different from 0. The other algorithms fail to decouple the linear contributions, since cross-entropy values different from zero can be seen in the off-diagonal elements for each one of the other methods, indicating the presence of some degree of coupling between input variable and the different estimated linear contributions. In addition, the larger cross-entropy values indicate stronger coupling between the variables. Only the outputs from the subspace system identification model produce larger crossentropy values between the input signals and their respective contributions when compared to the output from SIDE-ObSP. This might be attributed to the fact that the solutions provided by subspace system identification are less noisy. However, subspace system identification is not able to completely decouple the influence of other inputs onto each partial linear contribution.

### 4.2. Removal of Physiological Artifact from NIRS Signals

In this section we use SIDE-ObSP to decouple the influence of the variations in SaO<sup>2</sup> from the tissue oxygenation index (TOI), facilitating the use of TOI as a surrogate measurement of CBF for the assessment of cerebral Autoregulation (CA). TOI represents the ratio between oxygenated hemoglobin and total hemoglobin in the tissue. TOI is measured using spatially resolved spectroscopy as indicated in Suzuki et al. (1999). CA is a mechanism that tries to maintain a more or less stable CBF, despite the changes in MABP (Lassen, 1959). However, it has been shown that CA is not as simple as initially thought, and it comprises a set of responses from different mechanisms related to myogenic, neurogenic and metabolic actions (Panerai, 1998). Monitoring CA is important in order to avoid brain damage due


TABLE 1 | Cross-entropy values between each input variable, xi , and the estimated linear contributions, yˆ i , for ARX, ARMAX, adaptive filtering (Adap.), and subspace System identification (SS) models and SIDE-ObSP.

*The bold numbers indicate values of cross-entropy which are statistically different from zero.*

to ischemia and/or hemorrhage (Wong et al., 2008; Soul et al., 2007). TOI can be used as a non-invasive monitoring variable for CBF, it allows the continuous bedside monitoring of CA (Caicedo et al., 2011). However, TOI only reflects changes in CBF under a constant brain metabolism and constant arterial oxygen saturation (SaO2). In premature infants, even thought it is a strong statement, the first assumption can be considered valid during the first 3 days of life in the periods of analysis, which normally involve segments of 15 min. But, for the second assumption, the premature infants are likely to suffer from

TABLE 2 | Cross-entropy values between the input variables and the output generated by SIDE-ObSP.


TABLE 3 | Cross-entropy values between the input variables and the output generated by SIDE-ObSP for the complete studied population.


*The bold numbers indicate values of cross-entropy which are statistically different from zero.*

variations in systemic variables, especially in SaO2, during their stay at the neonatal intensive care unit (NICU). Under these conditions TOI cannot be used as a robust surrogate for CBF,

limiting the assessment of CA using NIRS in clinical practice. In this example we used data from 20 infants from the University Hospital Leuven (Belgium), with a gestational age of 28.4 ± 3.5 weeks and a birth weight of 1113 ± 499 g. Neonates were included following approval of the study protocol by the ethical board of the University Hospital, Leuven, Belgium, and after informed written consent was obtained from the parents. In all infants concomitant measurements of SaO2, measured by pulse oxymeter (Pulse oxymeter, Novametrix, USA), MABP, measured by an indwelling arterial catheter, and TOI, measured using spatially resolved NIRS (NIRO300, Hamamatsu,Japan), were obtained during the first 3 days of life with a sampling frequency of 1/3Hz. The total length of the recordings is between 6 and 9 h. The measurements were obtained during normal clinical care. For the signal decomposition algorithm, we used the measurements of SaO<sup>2</sup> and MABP as input variables, and the TOI as output variable. A selected segment of the recordings, where the influence of SaO<sup>2</sup> is clearly seen, is displayed in the first three panels in **Figure 3**. The figure shows a sudden drop in SaO<sup>2</sup> that is reflected in the TOI. We expect that SIDE-ObSP will produce as output the decomposition of TOI into 2 components, such that TOI <sup>ˆ</sup> <sup>=</sup> TOIMABP <sup>+</sup> TOISaO<sup>2</sup> , where TOISaO<sup>2</sup> is the component related to changes in SaO2, and TOIMABP is related to the changes in MABP, as shown in **Figure 4**.

The results from a representative segment are shown in the last panel of **Figure 3**. In order to evaluate that TOIMABP and TOISaO<sup>2</sup> are decoupled we used cross-entropy as explained in *Values given as median[min-max].*

the previous section. The cross-entropy values between MABP, SaO2, TOI, TOIMABP and TOISaO<sup>2</sup> for the representative subject are summarized in **Table 2**. As in the previous example, these values indicate that TOIMABP and TOISaO<sup>2</sup> are linked mostly to the dynamics of MABP and SaO2, respectively, and that there is not cross linked dynamics in the decomposed signals. According to these results we can consider that TOIMABP has been effectively decoupled from the variations in SaO2. In addition, based on the cross-entropy values, we can see that the TOISaO<sup>2</sup> has a stronger link to the changes in SaO2, than TOIMABP to the changes in MABP. This might indicate that the dynamic in SaO<sup>2</sup> affects more strongly the changes in TOI than the dynamic of MABP.

The results from the complete population are summarized in **Table 3**. It can be seen that cross-entropy values between SaO2/MABP and TOISaO<sup>2</sup> /TOIMABP are larger than the crossentropy values between SaO2/MABP and TOIMABP/TOISaO<sup>2</sup> . This indicate that the influence of other input variables in a given partial contribution are minimized. It is important to note that some of the cross-entropy values between an input and its respective partial contribution are low. Possible explanations for this can be that segments with a low variability in SaO<sup>2</sup> and MABP are not able to produce a proper description of the signals subspaces; that the segment under analysis is not coupled to the dynamics of the input signals, therefore the model fails to produce a component that is linearly related to the given input; or that the causal relationship, input vs. output, imposed in the model does not hold.

In the context of cerebral autoregulation assessment, TOIMABP can be used, directly, to assess the status of the CA

mechanism. This component not only represents a version of TOI decoupled from the variations in SaO2, but also represents the component of TOI that is linearly linked to the variations in MABP. SIDE-ObSP clearly offers the possibility of using NIRS for the assessment of cerebral autoregulation, even in the presence of changes in SaO2. This is an important result with a potentially high clinical impact that needs to be evaluated in further studies.

#### 4.3. Brain Hemodynamics Monitoring

In this section we use SIDE-ObSP in order to decompose the changes in cerebral and peripheral oxygenation into the partial contributions of each systemic variable, in order to evaluate the physiological responses caused by them, independently, in the peripheral and the cerebral circulation. We use a set of measurements obtained from a 3 years old infant under veno-arterial Extra corporeal Membrane Oxygenation (ECMO) procedure. ECMO is used to provide cardio-respiratory support to children with cardiac and/or respiratory problems. During this procedure the main vessels in the neck (right internal jugular vein and carotid artery) are cannulated. Blood is passed to an external oxygenator (via the vein cannula) and pumped back to the heart (via the arterial cannula). The heart and lungs are bypassed and allowed time to rest until they recover. Once there is indication that the patient's heart and lungs have recovered, they are being weaned off ECMO. During weaning, the ECMO blood flow is sequentially reduced, allowing the patients heart and lungs to take over. If during the weaning phase the patients can sustain normal function of their own heart and lungs they are decannulated and completely removed from ECMO. Patients under ECMO are at high risk of hemodynamic instability due to the possible alteration of the regulation mechanisms caused by multi-factorial reasons such as heparinitazion, hemodilution, and reduced arterial pulsatility (Annich et al., 2012).

These data was collected for research purposes, and the collection of the data was approved by the UCL, Institute of Child Health and Great Ormond Street Hospital for Children NHS Trust Research Ethics Committee. Written informed parental consent was obtained from the participants prior to inclusion. Cerebral tissue oxygenation (cTOI) was measured using a NIRS system (NIRO 200, Hamamatsu Photonics KK), using an optode located in the child's forehead; peripheral tissue oxygenation (pTOI) was measured in the calf using the same NIRS system. Concomitant measurements of MABP, end tidal CO<sup>2</sup> (EtCO2),

SIDE-ObSP. From top to bottom, partial contribution from HR, MABP, SaO2, EtCO2, and ECMO flow. On the right the frequency response of the individual subsystems linking each input with the output. This response was computed as the amplitude of the Fourier transform of the estimated impulse responses obtained from SIDE-ObSP. The frequency responses are organized presenting the one related to the cerebral hemodynamics first followed by the one related to the peripheral hemodynamics, for each input variable.

Heart Rate (HR), SaO2, and ECMO flow were also recorded during the ECMO weaning phase. Measurements were done during stepwise changes in the ECMO flow; the flow was reduced from baseline (100% ECMO flow) in steps of 10%, approximately every 10 min, until 70% of the baseline ECMO flow was reached, afterward the flow was increased back to baseline following the same profile (Caicedo et al., 2013b).

**Figure 5** displays some of the raw measurements taken during ECMO weaning. The upper plots in **Figure 5** shows the changes in HR, MABP, SaO2, EtCO2, and the ECMO Flow. The last plot presents the variations in cTOI and the pTOI.

The results are presented in **Figure 6**. In the left panel of the figure, the partial contributions of each systemic variable on the cerebral and peripheral circulation are shown. The right panel presents the frequency response for the different subsytems. By identifying the gain and the frequency band for each frequency response, the characteristics of the different mechanisms involved in the regulation of cerebral and peripheral hemodynamics can be compared. For instance, the magnitude of the frequency response can be used in order to estimate how much the changes on one of the systemic variables affect the corresponding hemodynamic variable. A larger magnitude in the frequency response, i.e., larger gain, indicates that the respective hemodynamic variables are more strongly affected by the specific systemic variation. This is also reflected in a partial contribution with a larger amplitude in the time domain. On the other hand, the bandwidth of the system can be computed by identifying the frequency region with a magnitude larger than 0, this information is useful in order to determine the dynamics of the partial contributions in time domain. Taking this information into account the results presented in **Figure 6** indicate that, in contrast with the peripheral circulation, the changes in cerebral hemodynamics caused by variations in MABP, HR, and ECMO flow are highly attenuated since the gain values are much smaller in the frequency response that correspond to the cerebral hemodynamics than the one from the peripheral hemodynamics. This is expected, since the regulation mechanisms that preserve brain hemodynamics are stronger than the regulation mechanisms in the peripheral hemodynamics. On the other hand, changes in SaO2, and EtCO<sup>2</sup> seem to affect cerebral and peripheral hemodynamics strongly, even though once again their impact on cerebral hemodynamics is more attenuated. This might indicate, to some degree, a similar regulatory action for changes in blood gases concentrations on both, the peripheral and the cerebral, vascular beds.

# 5. DISCUSSION

We have proposed a decomposition algorithm based on the use of ObSP, SIDE-ObSP, for the representation of NIRS signals into the partial linear contributions of different physiological signals. SIDE-ObSP assumes that the target signal can be modeled as the sum of the output of a set of mean average filters, one per each input variable. By approximating the structure of the model using MA filters, instead of an ARMAX structure, SIDE-ObSP is able to correctly define signal subspaces per each input variable, without contaminating them with the dynamics of other inputs. The decomposition is then achieved by properly defining oblique projectors that project the output signal onto the subspace spanned by a specific input variable. This oblique projector guarantees that the common dynamics between input variables is decoupled. Thereby, separating the partial linear contributions of each input, even if they are correlated.

SIDE-ObSP makes use of Tikhonov regularization using the derivative operator, in its matrix form, in order to produce smoother responses. The selection of the derivative operator was shown to be appropriate for the applications presented. However, other kernels can be easily integrated into the proposed algorithm, such as the tune-correlated (TC) kernel (Chen et al., 2011).

We demonstrated the performance of SIDE-ObSP with a synthetic example as well as in 2 medical applications in the field of cerebral hemodynamic monitoring. In the synthetic dataset we were able to retrieve the underlying impulse responses of the subsystems that generated the data, we compared the results obtained from SIDE-ObSP with other available methodologies and we showed that for this case, SIDE-ObSP outperformed them. In contracts with ARX, ARMAX, filter adaptive models, and subspace system identification, SIDE-ObSP was able to effectively decouple the dynamics between the input variables and the output, allowing to identify more accurately the underlying dynamics of the different subsystems. This is mainly due to the fact that the signal subspace of undesired input variables is contained in the Null subspace of the ObSP projector. In the first application example we use SIDE-ObSP in order to remove the influence of changes in SaO<sup>2</sup> from the changes in brain oxygenation measured using NIRS. This is a very important application since it facilitates the monitoring of cerebral autoregulation. In the second application example, we were able to provide a quantitative analysis of the differences between the influence of the systemic variables on the peripheral and cerebral hemodynamics thanks to the use of SIDE-ObSP.

Taking advantage of this decomposition scheme, we present the performance of SIDE-ObSP in 2 medical applications related to cerebral hemodynamics monitoring. First, we use SIDE-ObSP in order to filter out the physiological noise introduced by the variations of SaO<sup>2</sup> on the TOI. This is important for the non-invasive monitoring of cerebral autoregulation (Liem et al., 1995; Soul et al., 2007; Wong et al., 2008). We found that SIDE-ObSP was not only able to decouple the influence of the undesired fluctuations in SaO2, but it was also able to project the changes in TOI on the signal subspace spanned by MABP. This projection carries information about the linked dynamics between TOI and MABP, which can be directly used to quantify the status of the cerebral autoregulation mechanism. In addition, in the second application example, we showed that SIDE-ObSP decomposes the cTOI in the partial contributions of several systemic variables. We compared these results with the decomposition of pTOI and it was found, as expected, that the cerebral hemodynamics regulation mechanisms are able to mitigate and react appropriately to the changes in systemic variables in order to keep the brain hemodynamic homeostasis. The main advantage that SIDE-ObSP presents in this application field, is that it facilitates the individual interpretation and analysis of the influences of the changes in different systemic variables on the peripheral and cerebral hemodynamics. Such influences can be used for the analysis of the different mechanisms that are involved in the regulation of cerebral and peripheral hemodynamics. For instance, the frequency response of the sub-mechanism relating the changes in HR and pTOI acts like a band pass filter with a frequency response between 0 and 0.15 Hz, while the sub-mechanism relating changes in HR and cTOI acts as a low pass filter with a cut-off frequency around 0.1 Hz. This indicates that the partial contributions of HR on the cTOI in the time domain are smoother than the ones in the pTOI. However, response to changes in HR are 50 time stronger in the pTOI than in the cTOI, as indicated by the magnitude of the gain in the filters' pass band. Due to their frequency band we hypothesize this mechanism represents the sympathetic influence on the vascular tone of the cerebral capilar bed and the peripheral circulation. In a previous work we have studied the influence of the HR on the cerebral hemodynamics (Caicedo et al., 2014). Due to the assumption of a linear relationship between the HR and the cTOI and pTOI, no influence of the respiration and/or the cardiorespiratory coupling on the pTOI/cTOI has been included in this analysis, since the frequency range of the NIRS signals is restricted to a very low frequency range, we do not expect a large lineal influence from these variables. When comparing the frequency response of the subsystem linking MABP and TOI, it can be observed that the this subsystem behaves like a band-pass filter in the brain, whilst it exhibits a low-pass frequency response in the peripheral hemodynamics. Moreover, the response to changes in MABP are 20 times larger in the leg than in the brain. This indicates that changes in MABP are passively follow in the peripheral circulation, while only a small transient should be observed in the cerebral hemodynamic as response to changes in MABP, this transient behavior is typical from band pass systems, and it is caused by the fact that the DC values are highly attenuated, therefore any perturbation presented in the input should converge to zero. In addition, the frequency range of the frequency response in this subsystem indicates that the pTOI follows passively oscillations with a frequency smaller than 0.01 Hz (period larger than 100 s), while the frequency range that affects cTOI is restricted to frequencies between 0.01 and 0.1 Hz, representing oscillations between 10 and 100 s. Concerning the influence of SaO2, even though it affects both the cTOI and the pTOI, it can be seen that its influence is larger on the peripheral hemodynamics, represented by a gain 10 times larger in the frequency response. This response can be attributed to the larger compliance of the cerebral vascular bed, which is able to mitigate larger changes in CBF, and partially mitigate the effect of large changes in SaO<sup>2</sup> on the cTOI. The partial contributions of EtCO<sup>2</sup> to cTOI and pTOI are also similar; however, it can be seen that, in the brain, the changes in EtCO<sup>2</sup> are smoothed out stronger than in the peripheral circulation, this is also reflected in the lower cut-off frequency exhibit in the frequency response, around 0.01 Hz compared to 0.05 Hz respectively. This smoother behavior, like in the case of SaO2, can be caused by a higher capillary compliance in the brain than in the periphery. Finally, it is interesting to observe that the brain acts as a band-pass filter in the presence to changes in the ECMO flow, while the peripheral circulation exhibits a low-pass filter behavior. The protection mechanisms in the brain react to sudden changes in the ECMO flow which afterward are regulated, preserving a stable brain hemodynamics, whilst in the periphery these changes are reflected stronger and passively, due to the lowpass filter characteristic. In addition, the response to changes in ECMO flow in the pTOI are 80 time larger than the response in cTOI. These comparison and analysis was possible thanks to the use of SIDE-ObSP. Monitoring of the cerebral hemodynamics regulations mechanisms is critical for the prevention of brain damage and its consequent sequela, however, mathematical tools to properly monitor them are scarce, mainly, due to the complex linked dynamics of all the different mechanisms that are involved in cerebral hemodynamics regulation (Peng et al., 2008). In this context SIDE-ObSP represents one mathematical tool that can be effectively used in this field, in particular SIDE-ObSP can be useful in the study of the hemodynamic low frequency oscillations (LFO), these oscillations have been observed in since more than 150 years ago, but their origin and physiological explanation is still elusive. The complexity for the interpretation of the LFO is attributed to the coupling of different physiological processes (Sassaroli et al., 2012), it is here where SIDE-ObSP can be used in order to identify which systemic variations are more likely to be related to these oscillation, helping to unrevealed their possible origin and their potential link with cerebral autoregulation.

In addition, other algorithms can be used in order to decompose NIRS measurements in sources that might relate to physiological processes. In this context, Santosa et al. (2013) have developed an algorithm to remove noise sources, physiological or external, from functional NIRS (fNIRS) measurements using independent component analysis (ICA). ICA generally requires a set of measurements that are contaminated by common sources, fNIRS fits naturally within this framework, however, for single channel NIRS measurements this condition is not met. In such cases, single channel ICA can be used Davies and James (2007). However, this approach assumes that the sources are sufficiently disjoint in the frequency domain, which is a condition that cannot be imposed in physiological processes. Also the sources that are obtained using ICA might or might not be related to physiological mechanism. In comparison SIDE-ObSP is able to decompose single channels NIRS measurements and relate them to specific physiological measurements.

However, SIDE-ObSP presents some limitations. First of all, it assumes a linear relationship between the input variables and the output, which is likely not to be the case in biomedical systems. But, if the linear component has a strong contribution SIDE-ObSP can still produce relevant results, as proven with the application examples illustrated in the manuscript. Second, SIDE-ObSP is highly dependent on the available amount of information. This is due to the fact that the projectors for a given input variable are computed using information about the signal subspaces spanned by the other input variables. In case that one of these variables is not included, its subspace will be considered orthogonal to the input signal subspace, even if it is not. This might lead to unexplained dynamics on the decomposition residuals. In the context of clinical monitoring this might not be an issue, due to the availability of different biomedical measurements at the bed side. In addition, since the results are interpreted based on the input data available, residuals linked to other physiological sources can be carefully taken into account for further analysis and interpretation. Furthermore, by studying the dynamics of these residuals, other physiological variables that are linked to the dynamics of the output signal can be found. Finally, SIDE-ObSP assumes stationarity in the processed data segment. In the presence of nonstationarities, the computed model will produce a kind of "average" response. This problem can be mitigated by selecting data segments that are short enough to be considered stationary. However, the length of the processed segments should be large enough in order to guarantee that the number of rows in the regressor matrix **A** is larger that the number of columns. Approaches like subspace tracking methods can be used in order to process continuous stationary segments. This adaptation is useful for online monitoring systems, where nonstationarities are likely to occur but the changes in the system response between two consecutive segments are smooth.

### 6. CONCLUSION

We have presented a new algorithm, SIDE-ObSP, for the decomposition of NIRS signals into the linear partial contributions of a given set of input variables. SIDE-ObSP uses ObSP in order to effectively decoupled the dynamics from different signals. The application examples highlight its potential use in clinical practice, showing that it might help to enhance diagnosis, understanding of underlying (patho-)physiology, as well as to assess treatment. The main advantage of SIDE-ObSP is that it allows the individual interpretation of the relation between each input variable used in the model and the desired output variable. However, it is important to take into account that its performance relies on the quality and availability of the data used to train the model. Results should be interpreted with caution, taking into account the origin of the available measurements that have been used, and if these measurements are a direct reflection of the physiological variables of interest, or if they are surrogate variables that might include dynamics of different physiological systems. SIDE-ObSP can be extended to be used in other fields not exclusively for biomedical applications. SIDE-ObSP can be easily adapted for online monitoring by means of subspace tracking algorithms. Furthermore, due to its formulation, it can easily be adapted to a nonlinear regression framework based on kernels, such as LS-SVM regression models (Suykens et al., 2003).

## AUTHOR CONTRIBUTIONS

AC was involved in the conception and development of the algorithm presented in this manuscript. CV and BH provided intellectual support in the development of the methodology. MP and IT provided the data for the ECMO application study as well as intellectual support in the interpretation of the results. SV was in charge of the overall coordination of this study. All authors participated in the discussion of the results and the preparation, edition, and revision of the manuscript.

#### ACKNOWLEDGMENTS

This work is supported by: Bijzonder Onderzoeks- fonds KU Leuven (BOF): The effect of perinatal stress on the

#### REFERENCES


later outcome in preterm babies (# C24/15/036); Research Council KUL: GOA MaNet, PFV10002 (OPTEC), several Ph.D.postdoc and fellow grants; Flemish Government: FWO: Postdoc grants, projects: G.0427.10N (Integrated EEG-fMRI), G.0108.11 (Compressed Sensing) G.0869.12N (Tumor imaging) G.0A5513N (Deep brain stimulation). IWT: Ph.D. grants, projects: TBM 070713-Accelero, TBM 080658-MRI (EEG-fMRI), TBM 110697-NeoGuard. iMinds: SBO dotatie 2013, ICONs: NXT\_Sleep, FallRisk. Flanders Care: Demonstratieproject Tele-Rehab III (2012–2014). Belgian Federal Science Policy Office: IUAP P719 (DYSCO, "Dynamical systems, control and optimization," 2012–2017); ESA AO-PGPF-01, PRODEX (CardioControl) C4000103224. Belgiuan Foreign Affairs-Development Cooperation: VLIR UOS programs. EU: RECAP 209G within INTERREG IVB NWE programme, EU HIP Trial FP7-HEALTH 2007–2013 (n 260777), EU MC ITN TRANSACT 2012 (n◦ 16679), ERC Advanced Grant: BIOTENSORS (n◦ 339804), ERASMUS EQR: Community service engineer (n◦ 539642-LLP-1-2013). AC is a postdoctoral fellow of the Research Foundation-Flanders (FWO).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Caicedo, Varon, Hunyadi, Papademetriou, Tachtsidis and Van Huffel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Noninvasive Assessment of the Effect of Position and Exercise on Pulse Arrival to Peripheral Vascular Beds in Healthy Volunteers

Yurie Obata<sup>1</sup> , Qi J. Ong<sup>2</sup> , J. T. Magruder <sup>3</sup> , Helen Grichkevitch<sup>1</sup> , Dan E. Berkowitz <sup>1</sup> , Daniel Nyhan<sup>1</sup> , Jochen Steppan<sup>1</sup> and Viachaslau Barodka<sup>1</sup> \*

*<sup>1</sup> Division of Cardiac Anesthesia, Department of Anesthesiology and Critical Care Medicine, Johns Hopkins University School of Medicine, Baltimore, MD, USA, <sup>2</sup> Newcastle University School of Medicine, Newcastle, UK, <sup>3</sup> Division of Cardiac Surgery, Department of Surgery, Johns Hopkins University School of Medicine, Baltimore, MD, USA*

Background: The effects of position and exercise on pulse wave distribution across a healthy, compliant arterial tree are not fully understood. We studied the effects of exercise and position on the pattern of pulse arrival times (PATs) in healthy volunteers. Moreover, we compared the pulse arrival time ratios to the respective distance ratios between different locations.

Edited by: *Alexey Goltsov, Abertay University, UK*

#### Reviewed by:

*Ningjun Li, Virginia Commonwealth University, USA Kenji - Shigemi, University of Fukui, Japan*

> \*Correspondence: *Viachaslau Barodka vbarodk1@jhmi.edu*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *01 September 2016* Accepted: *17 January 2017* Published: *06 February 2017*

#### Citation:

*Obata Y, Ong QJ, Magruder JT, Grichkevitch H, Berkowitz DE, Nyhan D, Steppan J and Barodka V (2017) Noninvasive Assessment of the Effect of Position and Exercise on Pulse Arrival to Peripheral Vascular Beds in Healthy Volunteers. Front. Physiol. 8:47. doi: 10.3389/fphys.2017.00047* Methods: Thirteen young healthy volunteers were studied, using an electrocardiogram and plethysmograph to simultaneously record pulse wave arrival at the ear lobe, index finger and big toe. We compared the differences in PAT between each location at rest and post-exercise in the supine, sitting, and standing position. We also compared the PAT ratio (toe/ear, toe/finger, and finger/ear) to the corresponding pulse path distance ratios.

Results: PAT was shortest at the ear then finger and longest at the toe regardless of position or exercise status. PATs were shorter post-exercise compared to rest. When transitioning from a standing to sitting or supine position, PAT to the ear decreased, while the PAT to the toe increased, and PAT to the finger didn't significantly change. PAT ratios were significantly smaller than predicted by the path distance ratios regardless of position or exercise status.

Conclusions: Exercise makes PATs shorter. Standing position decrease PAT to the toe and increase to the ear. We conclude that PAT and PAT ratio represent the arterial vascular tree properties as surely as pulse transit time and pulse wave velocity.

Keywords: pulse arrival time, plethysmograph, body position, exercise, pulse wave velocity

# INTRODUCTION

When the heart contracts it ejects a bolus of blood, the stroke volume, into the arterial vascular system, which is then distributed to the peripheral tissues. The time period needed for a pulse wave to travel from the heart to the peripheral tissues depends on the distance traveled and the velocity of the pulse wave, which in turn is determined by various physical properties such as vessel diameter, wall thickness, and compliance (Bramwell and Hill, 1922).

**43**

Pulse wave velocity (PWV) has been studied extensively as a predictor for adverse cardiovascular events (Steppan et al., 2011). PWV is the ratio of vascular path length and the pulse transit time (PTT) between two sites of measurement. Pulse arrival to the peripheral vascular beds, the PTT, is the time from aortic valve opening till the foot of the pulse waveform at the peripheral site. Since it is impossible to directly measure aortic valve opening from an ECG alone, Pulse arrival time (PAT) has been widely used instead of PTT. However, these two are distinctly different. PAT is defined as the time delay between the peak of the R wave of the ECG waveform and the arrival (upstroke) of the arterial pulse wave in the periphery. Hence, PAT includes both pre-ejection time (ventricular electromechanical delay plus isovolumic contraction time) and PTT (Zhang et al., 2011).

Prior studies reported the PAT to different peripheral vascular beds in subjects of different age, height, and with different blood pressures (Allen and Murray, 2002; Nitzan et al., 2002). The PAT was consistently shortest at the ears, and longest at the toe. They showed that higher blood pressures and age were associated with shortening of the PAT, suggesting that PWV differs in distinctive peripheral vascular beds (Nitzan et al., 2002). Indeed Liu et al. showed that PWV is highest between the heart and toes, a slightly lower between the heart and fingers, and lowest between the heart and earlobe (Liu et al., 2011). In this manuscript, we studied the effects of position and exercise on the PAT to different parts of the body. Changes in body position triggers multiple cardiovascular responses due to the changes in hydrostatic pressure and sympathetic activity. In addition, exercise affects vascular properties through both vasodilation of arterioles and increased vasoconstriction due to sympathetic stimulation (Sharman et al., 2005; Roy and Secomb, 2014). Moreover, we compared the pulse arrival ratios to the respective distance ratios between different locations.

Specifically, we examined the PATs to three different vascular beds (ear lobe, index finger, and big toe) for the same heartbeat in three different positions (standing, sitting, and supine) at rest and post-exercise in young healthy subjects, as these factors might significantly affect PAT at different peripheral vascular beds.

### MATERIALS AND METHODS

This study was approved by The Johns Hopkins Medicine Institutional Review Board (IRB00074229). We enrolled 13 healthy volunteers with no history of vascular or cardiac disease, age 23–41 years old. Recruitment was done through e-mail or word of mouth in accordance with the Internal Review Board consent scripts. Inclusion criteria were: Healthy adults, age 21–50 years, both genders. Exclusion criteria were: Subject refusal to participate, known cardiovascular disease of any kind, pregnancy, and any disability preventing mild physical exertion. Two subjects who joined the study were excluded. The first subject was excluded due to the inability to finish the study protocol. The second subject was excluded as we were unable to obtain a plethysmograph signal on the toes. After verifying that each subject had no restrictions to participate in the study, we measured and recorded each subject's weight, half wingspan (distance from sternal notch to index finger with arm in 90◦ lateral extension), and self-reported height.

### Study Protocol

A standard three lead electrocardiogram (ECG; Bio Amp FE132, ADInstruments, Australia) was placed for continuous monitoring of electrical cardiac activity. Next, we placed capillary plethysmograph sensors [MLT1020EC IR Plethysmograph (ear), MLT1020PPG IR Plethysmograph (finger and toe), ADInstruments, Australia] on both left and right sides for each of the following locations: ear lobes in a standing position, index fingers in a sitting position with hands hanging free by their sides, and on the big toes lying in the supine position. We simultaneously recorded the ECG and plethysmograph bilaterally for each location (ear, index finger, big toe). Then, we simultaneously recorded the ECG along with the plethysmographs from one unilateral ear, finger, and big toe for 30 s each in the standing, sitting, and supine position. The ECG and plethysmograph sensors were then removed from the subjects and a blood pressure cuff applied to record blood pressure in the standing, sitting and prone position respectively. For the exercise part of the experiment the subjects were required to perform 30 squats. The ECG and plethysmograph sensors were reattached to the subject's left ear lobe, index finger, big toe, and the recording redone in the standing, sitting, and supine positions. All data for the "post-exercise" portion was collected within 3 min of the subject completing 30 squats. The ECG and plethysmograph signals were simultaneously digitally converted at 1 kHz (PowerLab 4/26, ADInsruments, Australia) and analyzed by eye and hand (LabChart 8.0, ADInstruments, Australia). The PowerLab has intrinsic filters built in the Bio Amp FE132 (ADInstrumens, Australia; see https://www.adinstruments.com/products/PowerLab). Since all signals were recorded simultaneously and through the same AD converter (PowerLab 4/26, ADInstruments, Australia) there should be no time delays according to the manufacture.

#### Data Extraction

From the data collected, the PAT to each location (ear lobe, index finger, and big toe) was assessed by calculating the time delay between the peak of the R wave on the ECG and the first subsequent positive inflection on the plethysmograph tracing. To compare PATs to different tissue beds from the same heartbeat, we took the corresponding R wave of the ECG, which was assigned a time value of zero. The representative data is shown in **Figure 1**. Periods consisting of 10 consecutive heart beats from each position (standing, sitting, supine) before and after exercise were then used to calculate the mean PAT for the three locations (ear lobe, index finger, big toe). To present the relative time differences between PAT to either ear, finger, or toe; we calculated the PAT ratios, toe/ear, toe/finger, and finger/ear. The collected data was then tabulated and used for statistical analysis.

**Abbreviations:** PAT, pulse arrival time; PWV, pulse wave velocity.

(second from the top), and ear (third from the top) and for 30 s each in standing, sitting, and supine positions. The PAT to each location (ear lobe, index finger, and big toe) was assessed by calculating the time delay between the peak of the R wave on the ECG and the first subsequent positive inflection on the plethysmograph trace. Periods consisting of 10 consecutive heart beats from each position (standing, sitting, supine) before and after exercise were then used to calculated the mean PAT for the three locations (ear lobe, index finger, big toe). ECG, electrocardiogram; PAT, pulse arrival time.

### Statistical Analyses

Data was analyzed using GraphPad Prism version 6.0 (GraphPad Software, San Diego, California, USA). Continuous variables are reported as mean ± standard deviation (SD). A paired ttest was used to compare the PATs on the left side with the PATs in the right side and to compare the PATs at rest with the PATs post-exercise at each location and in each position. Repeated-measures analysis of variance (ANOVA) with a Tukey test were used for multiple comparisons between PATs measured at different locations and positions as well as PAT ratios between different locations. Repeated-measures ANOVA with a Dunnett test were used for multiple comparisons between the path distance ratios and PAT ratios in different positions. The threshold for statistical significance was chosen to be P < 0.05.

### RESULTS

Demographics and baseline characteristics of the volunteers are presented in **Table 1**. Average age was 30 years and ranged from 23 to 41 years. Females made up 45% of the cohort and males 55%. At rest systolic blood pressure was 114 ± 13 mmHg and diastolic blood pressure 66 ± 7 mmHg.

PATs at each location (ears, toes, fingers) were not significantly different for the left and right side (p = 0.49 at ear, p = 0.22 at finger, p = 0.39 at toe; **Figure 2**). We therefore elected to compare only one side (left) to measure the PATs at different locations.

Mean and SD of the PATs and PAT ratios are tabulated in **Table 2**.

# Effects of Location and Exercise on Pulse Arrival (Figure 3)

In all positions, both at rest and post-exercise, the PATs between ear and toe were significantly different (shortest at the ear and longest at the toe).

Post-exercise PATs at all three locations and in all three positions were shorter compared to PATs at rest although some of them were not statistically significant.

# Effects of Position and Exercise on Pulse Arrival at the Ear (Figure 4A)

At rest, PAT to the ear was longest in the standing position (0.15 ± 0.02 s) compared to the sitting (0.14 ± 0.02 s, p = 0.03) and supine position (0.12 ± 0.02 s, p < 0.001). After exercise, PATs at the different locations were not significantly different (0.12 ± 0.01 s standing, 0.12 ± 0.01 s sitting, 0.11 ± 0.01 s supine).



*DBP, diastolic blood pressure; SBP, systolic blood pressure; SD, standard deviation.*

# Effects of Position and Exercise on Pulse Arrival at the Finger (Figure 4B)

The effect of position on pulse arrival to the finger was minimal compared to the ear or toe. Both at rest and post-exercise, the PATs measured in different positions were not significantly different (at rest: 0.21 ± 0.02 s standing, 0.20 ± 0.01 s sitting, 0.20 ± 0.01 s supine, post-exercise: 0.19 ± 0.03 s standing, 0.18 ± 0.01 s sitting, 0.18 ± 0.01 s supine).

# Effects of Position and Exercise on Pulse Arrival at the Toe (Figure 4C)

At rest, the PAT at the toe was longest in the supine position (0.29 ± 0.02 s) compared to standing (0.26 ± 0.02 s, p = 0.01) or sitting (0.26 ± 0.02 s, p = 0.03). The PATs in the standing and sitting positions were not significantly different (p = 0.98). Similarly, after exercise, the PAT was longest in the supine position (0.29 ± 0.02 s) compared to standing (0.25 ± 0.02 s, p < 0.001) or sitting (0.25 ± 0.01 s, p < 0.001; **Figure 4C**). The PATs in the standing and sitting positions were not significantly different (p = 0.46).

## Toe/Ear Pat Ratio and Toe/Ear Pulse Path Distance Ratio (Figure 5A)

Both at rest and post-exercise, the PAT ratio was significantly smaller than the pulse path distance ratio (5.7 ± 0.2), regardless of the position (p < 0.001).

At rest, the PAT ratio was smallest when standing (1.84 ± 0.20) compared to sitting (1.98 ± 0.27, p = 0.01) or supine (2.49 ± 0.32, p < 0.001).

After exercise, the PAT ratio was smaller when standing (2.10 ± 0.19) and sitting (2.11 ± 0.20) compared to being supine (2.55 ± 0.26, p < 0.001).

In the standing position, the PAT ratio was smaller at rest compared to post-exercise (p < 0.001). In the sitting and supine positions, the PAT ratios at rest were not significantly different compared to post-exercise (p = 0.09 and p = 0.49, respectively).

# Toe/Finger Pat Ratio and Toe/Finger Pulse Path Distance Ratio (Figure 5B)

Both at rest and post-exercise, the PAT ratio was significantly smaller than the pulse path distance ratio (1.7 ± 0.1), regardless of the position (p < 0.001).

At rest, the PAT ratio was smallest when standing (1.28 ± 0.10) compared to sitting (1.33 ± 0.11, p = 0.04) or supine (1.47 ± 0.32, p = 0.01).

After exercise, the PAT ratio was smaller when standing (1.32 ± 0.13) or sitting (1.38 ± 0.06) compared to supine (1.57 ± 0.16, p = 0.02: standing vs supine, p < 0.001: sitting vs. supine).

In both the standing and sitting position, the PAT ratios at rest were not significantly different compared to post-exercise (p = 0.33, p = 0.18, respectively). In the supine position, the PAT ratio was smaller at rest compared to post-exercise (p = 0.03).


*PAT, pulse arrival time.*

*P-values are bold when P* < *0.05.*

### Finger/Ear Pat Ratio and Finger/Ear Pulse Path Distance Ratio (Figure 5C)

Both at rest and post-exercise, the PAT ratio was significantly smaller than the pulse path distance ratio (3.3 ± 0.1), regardless of the position (p < 0.001).

At rest, the PAT ratio was smaller when standing (1.43 ± 0.11) and sitting (1.49 ± 0.19) compared to being supine (1.69 ± 0.21, p < 0.001).

After exercise, the PAT ratios for standing (1.60 ± 0.23), sitting (1.53 ± 0.15), and being supine (1.63 ± 0.15) were not significantly different (p = 0.37).

In the standing position, the PAT ratio was smaller at rest compared to post-exercise (p = 0.04). In the sitting and supine positions, the PAT ratios at rest were not significantly different compared to post-exercise (p = 0.29 in sitting, p = 0.31 in supine).

## DISCUSSION

The current study was performed in young healthy individuals to investigate the effects of position and exercise on the physiologic mechanisms of pulse wave distribution across a healthy, compliant arterial tree. Using noninvasive methods, we studied the pattern of pulse wave arrival at three distinct peripheral vascular beds at various positions and after exercise in healthy young volunteers.

Despite the fact that we used similar methodology as other groups, our study was focused on the effect of exercise and position on PAT, and we believe that our work is substantially different and adds new knowledge to the prior work published by Nitzan et al. (2002) and the group of Allen and Murray (2002). Moreover, to the best of our knowledge we are the first group that compared pulse arrival ratios to the respective distance ratios.

FIGURE 4 | The effect of position and exercise on the PAT. Mean and *SD* in PAT at rest and post-exercise in different positions. (A) Comparison of PAT measured at ear. (B) Comparison of PAT measured at finger. (C) Comparison of PAT measured at Toe. \*\*\**P* < 0.001; \*\**P* < 0.01; \**P* < 0.05; n.s., not significant. PAT, pulse arrival time; SD, standard deviation.

Similar to the studies mentioned above (Allen and Murray, 2002; Nitzan et al., 2002) we found that the pulse wave always arrives first at the ear, then at the index finger, and finally at the big toe—regardless of position or exercise status. This can be explained by the difference in the distances of the three locations from the heart. However, the PAT differences between distinct locations presented by the ratios of the PAT were relatively small compared to the distance ratios (**Table 1**). For example, the distance from the heart to the ear is on average 5.7 times shorter than to the toe, but the pulse arrives at the ear compared to the toe only 1.8–2.5 times sooner. This suggests that the arterial vascular tree actively maintains the smallest time difference in pulse wave arrival to the peripheral tissue beds irrespective of its distance from the heart. Moreover, the ratios are smaller for standing and sitting compared to the supine position, meaning that the PATs are more aligned in the vertical position compared to horizontal. We could speculate that a more simultaneous pulse arrival throughout the body is more important in the upright positions when people are likely to be more active, but we do not have any data to substantiate this claim.

Blood is transported along the arterial tree from the heart to peripheral tissues in order to supply those tissues' metabolic demand and to maintain homeostasis. This blood flow from the heart through the arterial system is provided in a pulsatile manner up to the pre-sphincter arterioles (Peterson, 1954). This allows the detection of a pulsatile waveform at the peripheral tissues, using plethysmography or pulse oximetry. Moreover, the time taken for the pulse wave to reach the vasculature of various peripheral tissue beds depends on a specific pulse propagation speed or pulse wave velocity (PWV) in each portion of the vascular tree (Bramwell and Hill, 1922). As one might infer, the distance of a particular vascular pathway always remains relatively constant, however PWV might undergo significant changes in particular arterial segments due to factors such as vascular tone and distending pressure (Bramwell and Hill, 1922). Since the arterial system does not contain valves, changes in position such as from lying to standing could hypothetically cause blood to pool in the lower extremities due to the formation of a hydrostatic gradient from the head to toe. In reality, in healthy individuals, the body produces a myogenic response, reacting to this change in hydrostatic pressure, by increasing the vascular tone through vasoconstriction in the lower extremities, leading to an increased distending pressure and wall tension in those arteries (Schubert and Mulvany, 1999). This increase in wall tension contributes toward an increased PWV and hence a shorter PAT (Bramwell and Hill, 1922; Learoyd and Taylor, 1966).

### Exercise

We hypothesized that the PAT shortens with exercise due to an increase in force of myocardial contraction, increased cardiac output, higher vascular wall tension, and central arterial wall stiffness (Strandell and Shepherd, 1967). Indeed our findings demonstrate that the PAT across all locations was shorter post-exercise. However, the pulse wave arrived relatively late at the lower extremity compared to the finger and ear as evidenced by increasing ratio of the toe/ear PAT and toe/finger PAT after exercise as compared to rest. One potential explanation is vasodilation in the metabolically active lower limbs and vasoconstriction in the less metabolically active upper limbs and head, to optimize metabolic/perfusion matching (Calbet and Joyner, 2010). Given that the subjects performed (squats) and therefore mainly utilized the muscles in the lower body, metabolic demand in the lower limb muscles would be higher than both the head and upper limbs. The net effect from the vasodilation in the lower limbs and vasoconstriction in the head and upper limbs was an increase in the toe/ear PAT ratios post-exercise. This is consistent with the fact that an increase in vascular tone (vasoconstriction) would increase PWV causing a decrease in PAT and vice versa (Bramwell and Hill, 1922; Sharman et al., 2005).

### Position

Our observations showed that a change in position from horizontal to vertical (supine to sitting or standing) led to a decrease in the PAT at the toe and an increase in the PAT at the ear. This might be explained by an increase in hydrostatic pressure in the arteries of the lower extremities and a decrease in the hydrostatic pressure in the arteries of the head and neck caused by the positional change from horizontal to vertical (Hasegawa and Rodbard, 1979). In healthy individuals, the vasculature adapts to the increase in hydrostatic pressure in the lower limbs by increasing wall tension or stiffness (Eiken et al., 2014). Without this response, blood would pool in the lower limbs and the perfusion of the upper limbs and head would be reduced. In many patients with orthostatic hypotension, disorders of the autonomic nervous system prevent this response from occurring, leading to reduced perfusion of the cerebral cortex and syncope (Stewart, 2002). Studies have shown that an increase in hydrostatic pressure and hence increased vascular wall tension causes an increase in PWV (Hasegawa and Rodbard, 1979). This explains our findings, wherein the PAT to the toes were shorter most likely due to the increase in wall tension and hence PWV, and the PAT to the ears were longer due to a decrease in wall tension and PWV when changing from a supine to standing position. However, overall, the PAT changes from lying supine to standing were not large. This might be due to fast adaptation by the vascular tone, autonomic nervous system signals, and myogenic responses as suggested previously (Pohl et al., 2000; Segal and Jacobs, 2001; Secomb, 2008; Roy and Secomb, 2014).

In our study, the pulse wave reached the ear before the index finger or big toe. This may be indicative of the body preferentially diverting blood to the cephalic region, although much more exploratory research is needed before coming to this conclusion. It is also unknown whether this pattern of pulse arrival is still true in individuals not included in the demographic studied (nonhealthy individuals, children, or the elderly). The data collected in this study could be used as a baseline to which future studies can be compared against.

Our results are consistent with published findings of others. Liu et al. used ECG and photoplethysmography to measure PWV between different sites of the body in supine position (Liu et al., 2011). They found that PWV is highest between the heart and toes, a little lower between the heart and fingers and lowest between the heart and earlobe. They also report that PWV of the left and right sides were the same, which implies that pulse wave arrival time is identical to each location between the left and right. Nitzan et al. at all found that pulse arrival time to the toe is longer compared to the finger in supine position, which also supports our observation (Nitzan et al., 2002). Work of Allen and Murray on healthy subjects of different age, height, and with different blood pressures showed that PATs become shorter with increasing age and systolic blood pressure (Allen and Murray, 2002). Similar to our results they reported that PAT is shorter at the ear compared to the finger and toe in supine position. Our work expanded prior knowledge by focusing on the effects of exercise and position on PATs. Moreover, we compared the PAT ratios to their respective distance ratios. Our results suggest that the arterial vascular tree properties minimize the difference in the PAT between different peripheral tissues despite large changes in path distance.

#### Limitations

The study was conducted on young and healthy volunteers and is therefore limited to this population. Our method of recruitment undoubtedly introduced a degree of sampling bias. The PowerLabs converter used was unable to distinguish between pressure generated by capillary filling and extremity movements; subjects had to stay very still while measurements were taken to ensure a readable trace. One subject who met all the inclusion criteria had to be excluded from the study because the equipment was unable to provide a readable toe plethysmograph trace. The PowerLabs hardware allowed us to record only four tracings at a time, such that we could not simultaneously record signals from both sides of the body. Rather, we had to do bilateral measurements first to confirm that each pulse arrived simultaneously to both sides at each level (toe, finger, and ear). For the exercise portion, we found it difficult to standardize the amount of physical activity based on our subjects various fitness and strength levels. Therefore, we felt squats would be an exercise that would be best suited to each individual's body weight and strength level. In the analysis of the data, we did not normalize for subjects' blood pressure or heart rate, which have been known to confound PWV (Lantelme et al., 2002). The choice of the peak of the R-wave on the ECG as the starting point of our PAT measurements also has its limitations as it includes the pre-ejection systolic phase of the ventricular contraction, the period between of isovolumetric ventricular contraction before the opening of the aortic valves (Geddes et al., 1981). However, the pre-ejection phase in healthy subjects is about 35 ms, whereas the PAT are in range of 150–300 ms, thus the potential error is in range of 15% (Biering-Sørensen et al., 2016). In future studies, the usage of a cardiac microphone to record the S1 heart sound could be more appropriate in obtaining a PAT sample.

### CONCLUSIONS

The main conclusions from this study are:

(1) The pulse wave always arrives at the ear first, then at the index finger, and finally the big toe—regardless of position or exercise status. This can be explained by the difference in the distances of the three locations from the heart.

# REFERENCES


Hence, we conclude that PAT and PAT ratio represent the arterial vascular tree properties as surely as pulse transit time and pulse wave velocity.

In the future, we anticipate studying PATs in aging and diseased populations, as well as with differential exercise regimens. It is known that both aging and cardiovascular disease significantly affect PWV, and the presenting mechanisms might be distorted in those groups. The current study could serve as the natural comparison group for future investigations in the elderly and in patients with cardiovascular disease.

#### AUTHOR CONTRIBUTIONS

Study design: VB. Data acquisition: QO, HG, VB. Data analysis: QO, YO, JM, VB. Interpretation of data, drafting the work, and revising it: QO, YO, DB, DN, JS, VB.


Zhang, G., Gao, M., Xu, D., Olivier, N. B., and Mukkamala, R. (2011). Pulse arrival time is not an adequate surrogate for pulse transit time as a marker of blood pressure. J. Appl. Physiol. 111, 1681–1686. doi: 10.1152/japplphysiol.0098 0.2011

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Obata, Ong, Magruder, Grichkevitch, Berkowitz, Nyhan, Steppan and Barodka. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Bifurcation in Blood Oscillatory Rhythms for Patients with Ischemic Stroke: A Small Scale Clinical Trial using Laser Doppler Flowmetry and Computational Modeling of Vasomotion

Alexey Goltsov <sup>1</sup> \*, Anastasia V. Anisimova<sup>2</sup> , Maria Zakharkina<sup>2</sup> , Alexander I. Krupatkin<sup>3</sup> , Viktor V. Sidorov <sup>4</sup> , Sergei G. Sokolovski <sup>5</sup> and Edik Rafailov <sup>5</sup> \*

#### Edited by:

Jincai Luo, Peking University, China

#### Reviewed by:

Jingyan Han, Boston University, USA Jun Li, Laboratoire Aimé Cotton (CNRS), France

#### \*Correspondence:

Alexey Goltsov a.goltsov@abertay.ac.uk Edik Rafailov e.rafailov@aston.ac.uk

#### Specialty section:

This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology

Received: 27 December 2016 Accepted: 02 March 2017 Published: 23 March 2017

#### Citation:

Goltsov A, Anisimova AV, Zakharkina M, Krupatkin AI, Sidorov VV, Sokolovski SG and Rafailov E (2017) Bifurcation in Blood Oscillatory Rhythms for Patients with Ischemic Stroke: A Small Scale Clinical Trial using Laser Doppler Flowmetry and Computational Modeling of Vasomotion. Front. Physiol. 8:160. doi: 10.3389/fphys.2017.00160 <sup>1</sup> Division of Science, School of Science, Engineering and Technology, Abertay University, Dundee, UK, <sup>2</sup> Department of Neurology, Neurosurgery and Medical Genetics, Pirogov Russian National Research Medical University, First City Hospital, Moscow, Russia, <sup>3</sup> Department of Functional Diagnostics, Priorov's Central Institute of Traumatology and Orthopedics, Moscow, Russia, <sup>4</sup> SPE "LAZMA" Ltd., Moscow, Russia, <sup>5</sup> Optoelectronics and Biomedical Photonics Group, Photonics and Nanoscience Group, Aston Institute of Photonic Technologies, Aston University, Birmingham, UK

We describe application of spectral analysis of laser Doppler flowmetry (LDF) signals to investigation of cerebrovascular haemodynamics in patients with post-acute ischemic stroke (AIS) and cerebrovascular insufficiency. LDF was performed from 3 to 7 days after the onset of AIS on forehead in the right and left supraorbital regions in patients. Analysis of LDF signals showed that perfusion in the microvasculature in AIS patients was lower than that in patients with cerebrovascular insufficiency. As a result of wavelet analysis of the LDF signals we obtained activation of the vasomotion in the frequency range of myogenic oscillation of 0.1 Hz and predominantly nutritive regime microcirculation after systemic thrombolytic therapy of the AIS patients. In case of significant stroke size, myogenic activity, and nutritive pattern microhaemodynamics were reduced, in some cases non-nutritive pattern and/or venular stasis was revealed. Wavelet analysis of the LDF signals also showed asymmetry in wavelet spectra of the LDF signals obtained in stroke-affected and unaffected hemispheres in the AIS patients. A mechanism underlying the observed asymmetry was analyzed by computational modeling of vasomotion developed in Arciero and Secomb (2012). We applied this model to describe relaxation oscillation of arteriole diameter which is forced by myogenic oscillation induced by synchronous calcium oscillation in vascular smooth muscle cells. Calculation showed that vasomotion frequency spectrum at the low-frequency range (0.01 Hz) is reciprocally modulated by myogenic oscillation (0.1 Hz) that correlates with experimental observation of inter-hemispheric variation in the LDF spectrum.

Keywords: laser Doppler flowmetry, acute ischemic stroke, cerebrovascular haemodynamics, blood flow oscillation, vasomotion modeling

# INTRODUCTION

Laser near-infrared spectroscopy (NIRS) and laser Doppler flowmetry (LDF) techniques along with spectral analysis of the optical signals become powerful method in clinical diagnostics of cerebrovascular diseases (Smith, 2011; Reinhard et al., 2003; Semyachkina-Glushkovskaya et al., 2016). High time resolution of NIRS and LDF methods allows measurement not only mean tissue oxygen saturation and blood flow (perfusion) but small hemocirculation fluctuations in cerebral oxygenation haemodynamic and metabolic changes (Mitra et al., 2016). These spectroscopic data were reported to bear valuable information on control of cerebral circulation in healthy and pathology (Schytz et al., 2010; Reinhard et al., 2014). One of directions in the analysis of medical optical spectra aims at the investigation of cerebrovascular autoregulation (CA) and its impairment as a result of cerebrovascular diseases. Using of NIRS and LDF methods in clinics allows non-invasive measurement of CA instead of invasive common techniques (Fantini et al., 2016) and the development of the multifunctional, inexpensive portable technologies enabling clinicians to easily perform continuous monitoring of patients during diagnosis, treatment, and posttreatment periods (Rafailov et al., 2015).

NIRS and LDF spectra contain information on temporal variation of the mean arteria blood pressure (MABP), cerebral blood flow (CBF), and oxygenation. They include both passive components such as heart and breathing pulses and the active ones such as spontaneous oscillation in low-frequency oscillation range from 0.05 to 0.15 Hz (LFO) and very low low-frequency oscillation ranged up to 0.05 Hz (VLFO) bands (Schytz et al., 2010). Specific frequency bands in optical spectra were reported to link to CA which is a critical mechanism to maintain stable blood flow in an ischemic liaison area. A change in CA is considered as an important factor in the diagnosis and treatment of cerebrovascular diseases (Reinhard et al., 2012; Liu et al., 2015). The most frequently used parameters for analysis of CA are phase shift (PS) in the frequency domains of LFO or VLFO and correlation coefficients in the time domains between MABP/CBF and oxygenated (oxyHb)/deoxygenated (deoxyHb) hemoglobin (Schytz et al., 2010; Reinhard et al., 2014). Clinical investigation of ischemic patients showed an increase in the correlation coefficients obtained for ischemic affected and nonaffected hemispheres that indicates weakening CA (Reinhard et al., 2012).

To extract information on cerebrovascular dynamics and CA abnormalities from NIRS and LDF signals, different techniques of spectral analysis are used, among which are Fourier transformation (Rossi et al., 2007; Vermeij et al., 2014), wavelet analysis (Stefanovska, 1999; Addison, 2015), signal correlation analysis (Toronov et al., 2013), and signal linear and nonlinear decomposition methods (Iatsenko et al., 2015; Caicedo et al., 2016). For further development and translation of the spectral analysis of NIRS and LDF signals as a clinical method, it will need detailed investigation of the molecular and physiological mechanisms behind LFO and VLFO in healthy and pathology (Schytz et al., 2010). Although spectral methods are intensively used in investigation of CA in pathological conditions, the link of these active oscillations with cerebrovascular diseases is not sufficiently explored. At present, investigation of a link between LFO amplitude and cerebrovascular diseases showed a significant response in spectral characteristics to such pathological conditions as carotid stenosis (Hu et al., 1999; Schytz et al., 2010), acute ischemic stroke (AIS) (Phillip and Schytz, 2014), traumatic brain injury (Chernomordik et al., 2016), and cerebrovascular occlusive disease (Phillip et al., 2013).

Intensive investigation of the spontaneous oscillation in microvascular flow, referred as vasomotion (oscillation of vascular diameter), revealed that its spectral characteristics are intrinsic hallmarks of microvascular regulation of blood flow in different tissues and organs (Kvandal et al., 2003). As a result of the spectral analysis of NIRS and LDF signals, several inherent frequency regions were identified and attributed to different physiological mechanisms of vascular control: (I) 0.6–2 Hz related to cardiac activity, (II) 0.145–0.6 Hz corresponds to respiratory activity, (III) 0.052–0.145 Hz linked to microvessel smooth muscle cell activity, (IV) 0.021–0.052 Hz connected to microvessel innervation, and region (V) 0.005 Hz– 0.0095 Hz associated to endothelia activity (Stefanovska, 1999; Rossi et al., 2008). Experimental investigation of vasomotion are accompanied by detailed computational modeling of its molecular mechanisms (Ursino et al., 1996; Koenigsberger et al., 2008; Arciero and Secomb, 2012). Results of the experimental and theoretical investigation are assumed that different mechanisms modulate vascular smooth muscle (VSM) activity in the corresponding frequency regions of the vasomotion spectrum. Detailed investigation of the vasomotion spectrum showed a significant response of the spectral components to pathophysiological conditions such as malignant melanoma (Lancaster et al., 2015), surgical trauma and anesthesia (Schmidt-Lucke et al., 2002), arterial hypertension, skin post-ischaemic reactive hyperaemia (Rossi et al., 2007, 2008), and experimental intervention such as the endothelium-dependant vasodilators (Kvandal et al., 2003). Clinical studies demonstrated that the spectral analysis of LDF signals in the specific frequency regions is a promising method of monitoring microvessel regulation in patients with cerebrovascular diseases (Rossi et al., 2008; Krupatkin and Sidorov, 2016).

In this work we investigated a feasibility of spectral analysis of LDF signals to evaluate a function status of cerebral microcirculation system in patients with acute and chronic disorder of cerebral circulation (Krupatkin, 2005). The developed method is based on the wavelet spectral analysis of an alteration in spectral components within different frequency regions of LDF spectrum discussed above. We also applied a computational model of vasomotion developed in Ursino et al. (1996) and Arciero and Secomb (2012) to analyse possible mechanisms underlying a change in LDF spectra observed in AIS patients.

### METHODS

Patient investigation was carried out in the Regional Vascular Centre of the Pirogov Russian National Research Medical University (RNRMU), Moscow, Russia. Clinical trial included 19 patients, 10 man and 9 women, aged from 27 to 76 years with the ischemic monohemispheric stroke and 12 of them were treated with thrombolytic therapy. The trial was carried out in accordance with the Declaration of Helsinki principles and approved by the Ethics Committee of the Pirogov Russian National Research Medical University. All patients participated in this trial have given the full consent on measurements in written and been informed of their right to discontinue participation at any time.

An initial average National Institute of Health Stroke Scale (NIHSS) score of the patients was 10.2 ± 4.7. Sixteen patients were classified with large-vessel disease and 13 patients with concurrent or undetermined causes according to TOAST criteria (Adams et al., 1993). Exclusion criteria in this trial were hemorrhagic stroke, cancer, infectious diseases including infections of the nervous system, mental disease, autoimmune disease, drug addiction, chronic alcoholism. Patients with AIS were evaluated by course of neurological, psychoemotional state, cognitive function, and laboratory results. The main characteristics of the patients are shown in **Tables 1**, **2**. As a reference cohort, 22 patients with chronic cerebrovascular encephalopathy were evaluated. Nineteenth percentage of patients in this group had a combination of arterial hypertension with atherosclerosis and in the remaining group of patients, arterial hypertension was more likely to be a leading factor because it was developed in young age.

LDF investigation was performed from the third to seventh day after the onset of ischemic stroke. Measurements were carried out on forehead in the left and right supraorbital regions over brow skin areas corresponding to a branch of the internal carotid (ophthalmic) arteries, which supplies blood to the supraorbital artery from internal carotid artery. LDF measurement was carried out using multifunctional laser device LAKK-M (SPE LAZMA Ltd., http://www.lazma.ru; Dunaev et al., 2014; Rafailov et al., 2015). LDF signals were measured up to 300 s with time step 1t = 0.05 s. Cerebral microcirculation parameter was obtained as a value of the LDF signals in arbitrary perfusion units (pu). Then LDF signals were analyzed using wavelet transformation to investigate spectral components in the specific frequency regions of cerebral microcirculation oscillations (Kvernmo et al., 1998). Statistical analysis included calculation of the means and standard deviations of the values of arbitrary perfusion and amplitudes of the main peaks in the wavelet spectra.

### Analysis of Wavelet Spectrum of LDF Signal

Spectral analysis of the LDF signal was carried out by the wavelet method with further determination of leading components in the specific frequency ranges of the spectrum. Following the method (Kvernmo et al., 1998), continuous wavelet transform of the LDF signal S(t) was calculated as

$$W(s, t\_0) = \frac{1}{\sqrt{s}} \int\_{-\infty}^{\infty} \mathcal{S}(t) \,\psi\left(\frac{t - t\_0}{s}\right) dt\tag{1}$$

where ψ t−t<sup>0</sup> s is the mother wavelet function which was chosen in the Morlet wavelet form

$$\psi\ (\nu) = \pi^{-1/4} \exp\left(-i2\pi f\_0 \nu\right) \exp(\nu^2/f\_1) \tag{2}$$

where f<sup>0</sup> and f<sup>1</sup> are dimensionless parameters, defining the wavelet center frequency and wavelet bandwidth respectively and v = t−t<sup>0</sup> s is a dimensionless variable with the wavelet scales s and t0. Finally, we calculated the wavelet power spectrum |W (s, t0)| 2 .

#### Fourier Analysis of LDF Spectrum

Spectral characteristics of the LDF signals were also analyzed using spectral density of auto-correlation function of the LDF signals which was calculated as

$$R\left(\tau\right) = \frac{\sum \left(\mathbf{S}(t) - \mathbf{\bar{S}}\right) \left(\mathbf{S}(t-\tau) - \mathbf{\bar{S}}\right)}{R\left(0\right)},\tag{3}$$

where S¯ is the average LDF signal over the measurement time T. The power spectral density (PSD) was calculated as discrete Fourier transformation of the auto-correlation function of the LDF signal R (τ ):

$$R(f\_k) = \sum R(\tau)e^{-i2\pi k \frac{n}{N}},\tag{4}$$

which is defined for discrete frequencies f<sup>k</sup> = k T , where k = N − 1. PSD was calculated as

$$P\left(\left\{k\right\} = \frac{1}{N^2} \left| R\left(\left\{k\right\}\right|\right|^2. \tag{5}$$

#### Computational Modeling of Vasomotion

To analyse the mechanism underlying the observed changes in the LDF spectra measured in AIS patients, we applied a nonlinear model of spontaneous (vasomotion) oscillation of microvessels which was developed in a series of theoretical papers (Ursino et al., 1996; Koenigsberger et al., 2008; Arciero and Secomb, 2012). The developed phenomenological models describe oscillation of arteriole diameter and blood flow under the vessel wall tension and share stresses. In these models, vasomotion involves oscillating contraction of the vascular smooth muscles (VSMs) in the arteriole wall which results from interaction between the mechanical tension of the vessel wall and the dynamics of tone generation in VSMs. In this work, we used a version of the vasomotion model developed in Arciero and Secomb (2012) and generalized it to consider myogenic oscillations of the arteriole wall which is suggested to be induced by synchronous calcium oscillation of VSMs. A kinetic model of this molecular mechanism of myogenic oscillations was developed in Koenigsberger et al. (2005). In our model, myogenic oscillations were phenomenologically considered as forcing tension induced by VSMs applied to the microvessel wall.

The model of spontaneous vasomotion of arteriole diameter (Arciero and Secomb, 2012) includes two ordinary differential equations (ODEs) which describe: (1) a change in diameter D of



AH, arterial hypertension; AS, atherosclerosis; CE, cerebral embolism; AT, arterial thrombosis.

TABLE 2 | Severity of patient status with hemisphere insult according to NIHSS score.


TABLE 3 | Parameters of the vasomotion model taken from Arciero and Secomb (2012).


an arteriole and (2) the response of vascular smooth muscles A (vascular tone) to the tension in the wall T and wall share stress τ:

$$\begin{cases} \frac{dD}{dt} = \frac{1}{\tau\_d} \frac{D\_0}{T\_0} \left( T(D) - T\_{tot}(D) \right) \\\frac{dA}{dt} = \frac{1}{\tau\_d} \left( A(D) - A\_{tot}(D) \right) \end{cases} \tag{6}$$

Here T (D) = PD/2 is the wall tension due to intravascular pressure P according to Laplace law. Ttot (D) is the tension generated by VSMs in the vessel wall which is defined by its passive Tpass(D) and active Tact(D) components (Arciero and Secomb, 2012):

$$T\_{tot}(D) = T\_{\text{pass}}(D) + A(D)T\_{act}(D) \,. \tag{7}$$

where

$$T\_{\text{pass}}(D) = C\_{\text{pass}} \exp\left(C\_{\text{pass}}^{'}\left(\frac{D}{D\_0} - 1\right)\right),\tag{8}$$

$$T\_{act}(D) = C\_{act} \exp\left(-\frac{D/D\_0 - C\_{act}'}{C\_{act}'}\right). \tag{9}$$

Activation function A(D) varies between 0 and 1 and takes a steady-state value Atot which is defined as a sigmoidal function of a stimulus S:

$$A\_{tot}\text{ (S)} = \left(1 + e^{-S}\right)^{-1},\tag{10}$$

where

$$\mathcal{S} = \mathcal{C}\_{myo}\mathcal{T} - \mathcal{C}\_{shear}\mathcal{T} + \mathcal{C}'\_{tone}.\tag{11}$$

The first term in function S corresponds to the myogenic response to the wall tension. The second term describes the sheardependence response and Ctone is a constant. Tension of VSMs T is produced by intravascular pressure and the wall share stress τ results from blood flow in small vessels. Constants D<sup>0</sup> and T<sup>0</sup> are normalized values of diameter and tension which were taken equal to passive reference vessel diameter and tension value T at D = D<sup>0</sup> (**Table 3**). The time scales τ<sup>d</sup> and τ<sup>a</sup> were selected to reproduce satisfactory low-frequency spontaneous oscillations in the range of 0.01 Hz observed in the LDF spectrum. The obtained values τ<sup>d</sup> = 5 s and τ<sup>a</sup> = 60 s are close to those given in Arciero and Secomb (2012).

We supposed that myogenic oscillations of small arteries and arterioles containing VSMs are generated by synchronous calcium oscillations in smooth muscle cells (SMCs) which force oscillations of vessel diameter (Koenigsberger et al., 2008). To describe these myogenic oscillations of the vascular diameter, we introduced additional tension T<sup>m</sup> into Equation (8) in the form of forcing oscillatory function of amplitude A<sup>m</sup> and frequency fm:

$$T\_m = A\_m \sin\left(2\pi f\_m t\right). \tag{12}$$

Frequency of the myogenic oscillation in the model was taken f<sup>m</sup> = 0.15 Hz that corresponds to the frequency region where myogenic oscillation was observed (Stefanovska, 1999).

In the model, we suggested that oscillations of the LDF signals result from the oscillations of the platelet velocity induced by active oscillations of the arteriole radius. Arteriole radius oscillations cause fluctuations in the hydraulic viscous resistance of an arteriole which are inversely proportional to the fourth power of the arteriole radius r according to Hagen-Poiseuille law, R = 8ηL/r 4 , where L is the length of a vessel segment and η is blood viscosity. The oscillations of the hydraulic viscous resistance cause oscillations in blood flow in the vessel segment which is defined by Poiseuille law,

$$Q = \frac{\pi \,\Delta P}{R},\tag{13}$$

where 1P is the pressure difference across the arteriole segment. Thus, oscillation of the LDF signals, which is proportional to Q, is defined by the active oscillation of the arteriole radius r. In this work, we compare spectra of the LDF signals of blood flow measured in patients and spectra of the arteriole radius

Frontiers in Physiology | www.frontiersin.org

oscillations calculated in the model at the different physiological conditions.

#### RESULTS

LDF investigation was carried out in patients with stroke localization in the right (10 patients) and left (9 patients) cerebral hemispheres. According to LDF data, the patients after thrombolytic therapy had asymmetric haemodynamic responses measured on the right and left sides of forehead corresponding to affected and unaffected hemispheres. Two typical LDF signals obtained on both forehead sides in the ischemic patients with a small focal lesion is shown in **Figures 1A,B**. The LDF signals on affected and unaffected hemispheres differ between each other both in amplitude and frequency. Comparison of the LDF signals in patients with different pathology showed significantly lower perfusion value (p < 0.05) in AIS patients (11.7 ± 5.5 pu and 13.3 ± 7.4 pu in affected and unaffected hemispheres respectively) than that in patients with chronic cerebrovascular encephalopathy (26.2 ± 4.4 pu).

Wavelet analysis was performed for five frequency ranges typical for LDF spectra of vasomotion: endothelial (0.0095– 0.02 Hz); neurogenic (0.02–0.06 Hz); myogenic (0.06–0.16 Hz); breathing (0.16–0.4 Hz) and heart rate (0.4–1.6 Hz). The results of the wavelet analysis of the LDF signal for the patient 1 with a small focal lesion (**Figure 1A**) are presented in **Figures 2A,B** for affected and non-affected hemispheres, respectively. Amplitude and normalized amplitude of LDF wavelet spectrum corresponding to the heart rate and myogenic frequency bands were A<sup>p</sup> = 0.36 and Am/σ = 0.57 for affected hemisphere and A<sup>p</sup> = 0.36 and Am/σ = 0.44 for unaffected one respectively, were σ is the standard deviation of the amplitude of LDF spectrum in this region. To estimate differences in the spectral characteristic of the wavelet spectra for affected and unaffected hemispheres, we calculated the means and spreading of the normalized wavelet spectral components (A/σ) and showed this comparison in Figure S1 in Supplement.

LDF wavelet spectra measured on a side of affected hemisphere in the AIS patients after thrombotic therapy showed that the heart rate rhythm has dominated among the passive oscillations in the spectrum that indicates blood inflow to microvessel bed. Analysing the breathing rhythm frequency spectrum, we did not observe significantly high amplitude values that might occur at venular stagnation of blood and a decrease in perfusion pressure. The amplitude of the breathing rhythm did not exceed the amplitude of the heart rate rhythm of 83.3% of patients. Typical feature of the LDF wavelet spectra observed in 83.3% patients was a considerable contribution of myogenic oscillation in the frequency range of active oscillations. This can indicates the myogenic component in autoregulation of the cerebral haemocirculation and suggests pre-capillary activity and an increasing number of open capillaries (Krupatkin, 2005). Thus, enhancement of myogenic activity of arteriole smooth muscles and boost of nutritive blood flow occurred after thrombolytic therapy (Krupatkin, 2005). We have to note, that this feature of LDF wavelet spectrum was observed in only 28.5% of the patients with ischemic stroke who did not receive thrombolytic therapy.

The analysis of the LDF signals measured on a side of affected hemisphere in patients with small (<5,000 mm<sup>3</sup> , group 1, n = 5) and with large (>5,000 mm<sup>3</sup> , group 2, n = 14) focal brain lesion mapped by CT was carried out. No significant difference in the perfusion values between groups 1 and 2 was observed, while distinct difference in wavelet spectral characteristics i.e., myogenic rhythm was pronounced for patients of the group 1. This clearly indicates activation of myogenic type of the autoregulation in the resistive arterioles and the nutritive mode of microcirculation (Krupatkin, 2005). The other active rhythms in the LDF wavelet spectrum were not well distinguishable. Amplitude ratio of heart rate to breathing rhythm was Ap/A<sup>b</sup> >1. For patients of the group 2, the other active oscillatory rhythms apart from myogenic oscillations were found, namely neurogenic and rare endothelial ones, that suggests shunt component of circulation and less pronounced nutritive mode of the cerebral blood microcirculation appeared to be active (Krupatkin, 2005). In some cases, myogenic oscillations in the spectrum were not dominated. In the passive oscillation range of the spectrum, the ration of the amplitudes was Ap/Ab<1 that testifies to venular stagnation (Krupatkin, 2009).

For further evaluation of the spectral characteristics of the LDF signals we analyzed spectral density of the auto-correlation function of the LDF signals [Equation (3) in Methods]. Autocorrelation functions R(τ ) of the LDF signals measured on unaffected (red line in **Figures 3A,B**) and stroke-affected (black line in **Figures 3A,B**) branches of the supraorbital arteries in two patients (**Figures 1A,B**) are given in **Figures 3A,B**. As seen, the auto-correlation functions of the LDF signals corresponding to stroke-affected and unaffected hemispheres differ from one another in their shapes and these differences are similar for both ischemic patients with a small focal lesion. Power spectral density of the auto-correlation function of LDF signals P(f<sup>k</sup> ) [Equation (5)] in the range of 0.01–0.1 Hz and in the heart rate range of 1 Hz are shown in **Figures 4A,B** respectively for the first patient (corresponding LDF signals and auto-correlation functions shown in **Figures 1A**, **3A**). Both auto-correlation function R(τ ) and its power spectral density P(f<sup>k</sup> ) calculated for unaffected (red lines) and stroke-affected (black lines) hemispheres significantly differ between each other (**Figures 3**, **4**). Comparison of these results with that of the LDF wavelet analysis showed similar differences of the LDF spectra measured on unaffected and AISaffected hemispheres.

To analyse the mechanism underlying the observed asymmetry in the LDF spectra in the left and right hemispheres in patients with mono-hemispheric stroke, we applied a nonlinear model of spontaneous (vasomotion) oscillation of microvessels which was developed by Arciero and Secomb (Arciero and Secomb, 2012; see Methods). This phenomenological model describes spontaneous oscillation of the arteriole diameter and blood flow under the wall tension and wall share stresses (Arciero and Secomb, 2012). In this work, authors defined a set of model parameters at which spontaneous (limit-cycle) oscillation of vascular diameter takes place and carried out bifurcation analysis of this oscillatory system (**Table 3**). Taking

FIGURE 1 | LDF signals measured on a side of stroke affected (black line) and unaffected (red line) hemispheres of patients 1 (A) and 2 (B) with a small focal AIS lesion.

small focal AIS lesion (Figure 1A). Five physiological frequency regions are marked: (I) cardiac activity, 0.6–2 Hz; (II) respiratory activity, 0.145–0.6 Hz; (III) myogenic activity of VSMs, 0.052–0.145 Hz; (IV) microvessel innervation, 0.021–0.052 Hz; (V) endothelial activity, 0.0095–0.021 Hz.

the same set of the model parameters (**Table 3**), we obtained spontaneous oscillations of vascular diameter in the frequency range of 0.01 Hz in this model (see dot line in **Figure 5A**).

Application of this model to the description of arteriole diameter oscillations in the presence of myogenic oscillations showed that the character of vasomotions significantly depends on the amplitude of myogenic oscillations. **Figure 5** shows results of the calculation of arteriole diameter oscillation, autocorrelation function and its power spectrum. Calculations were carried out at low amplitude A<sup>m</sup> = 1 kdyn/cm (**Figures 5A–C**) and high amplitude A<sup>m</sup> = 1.8 kdyn/cm (**Figures 5D–F**) of the myogenic oscillations. At low amplitude Am, the diameter oscillations were obtained to be quasiperiodic, i.e., high frequency oscillations were modulated by low frequency

spontaneous oscillations of 0.01 Hz (**Figure 5A**). At high amplitude of myogenic oscillations Am, high frequency oscillations prevail in the spectrum with less pronounced modulation by the low frequency of 0.02 Hz (**Figure 5D**). At higher amplitudes A<sup>m</sup> more than 2 kdyn/cm, the modulation by high frequency disappeared because of locked modes in the forced relaxation oscillation (data not shown).

Comparison of the spectrum obtained for low and high amplitudes of myogenic oscillation revealed that amplitude of low frequency limit-cycle oscillations (0.01 Hz) decreased when the amplitude of myogenic oscillation increased (see **Figure 5**). These changes in spectral density is in a satisfactory agreement with experimental data on asymmetry obtained for LDF spectra, i.e., a reduction of the myogenic rhythm component in the frequency range of 0.1 Hz is accompanied by a rise of the endothelia rhythm component in the region of 0.01 Hz in the stroke-affected hemisphere compared with the unaffected one (see **Figure 4**). Comparison of the experimental and theoretical auto-correlation functions showed that the model is capable to satisfactorily describe the change in their shapes observed at the increasing myogenic rhythm component that occurred in the stroke-affected hemisphere (**Figures 3**, **5B,E**).

To investigate the mechanism underlying the change in the spectrum of vessel diameter oscillations, we have calculated phase portrait of the vasomotion model [Equation (6)] and presented it at a diameter-activation (D/A) plane along with the oscillatory trajectory of the solution (**Figure 6**). In the absence of the myogenic oscillations, the trajectory of the solution corresponds to the limit-cycle oscillations (**Figure 6A**). These oscillations of vessel diameter at low frequency (0.01 Hz) has a specific sawshape character that corresponds to the 2D limit-cycle at the phase portrait (see the corresponding solution in **Figure 5A**, dotted line). At low amplitude of the myogenic oscillations the limit-cycle oscillations of the vessel diameter occur at low frequency of 0.01 Hz with superimposed high-frequency oscillations of 0.1 Hz (**Figure 6B**). At an increased amplitude of the SMC oscillations, the blowout of the 2D oscillatory trajectory from the stable branch of D-nullcline occurs (**Figures 6C,D**). At a further increase in amplitude of the SMC oscillation, the 2D trajectory of the limit-cycle oscillations degenerates into a semi-one dimensional trajectory that represents high-frequency oscillations slightly modulated by low-frequency oscillations (**Figure 6D**).

Summarising these results, we concluded that the changes in oscillation spectrum obtained in the model results from interaction of high-frequency myogenic oscillations and lowfrequency limit-cycle oscillations which are controlled by vascular tone. Increasing amplitude of myogenic oscillation causes "distortion" of a limit-cycle, weakening of low-frequency components in oscillation and prevalence of high-frequency myogenic oscillation in the vasomotion spectrum. This leads to subharmonic entrainment and locked mode in the forced relaxation oscillators. These phenomena were comprehensively investigated for the forced van der Pol relaxation oscillator (Jordan and Smith, 2007).

#### DISCUSSION

LDF investigation of patients with AIS was performed from 3 to 7 days after the onset of ischemic stroke on forehead in the supraorbital region, which supplies blood to the supraorbital artery from internal carotid artery. Analysis of the LDF signals showed that the blood perfusion in the microvasculature in the projection of the affected vascular field in AIS patients was lower compared to patients with cerebrovascular insufficiency. Wavelet analysis of LDF signals of the AIS patients revealed the blood flow increase in microcirculatory and a raise in myogenic activity of muscle-containing arterioles that predominantly indicate nutritive regime of the microcirculation after systemic thrombolytic therapy. In case of significant stroke size, the reduction of the myogenic activity and nutritive pattern of the microhaemodynamics were observed, and in some cases nonnutritive pattern and/or venular stasis were detected.

Results of the wavelet analysis of LDF signals also showed asymmetry in the LDF spectra detected in stroke-affected and unaffected hemispheres. A mechanism underlying this asymmetry was analyzed by the theoretical model of microvessel wall oscillations (vasomotions) developed in a series of publications (Ursino et al., 1996; Koenigsberger et al., 2008; Arciero and Secomb, 2012). Here we have expanded the model suggested by Arciero and Secomb (Arciero and Secomb, 2012) and applied it to describe interaction between myogenic oscillation, which as it was suggested to be induced by synchronous calcium oscillations in SMCs, and spontaneous oscillations coming from coupling regulation between endothelia cells (ECs) and SMCs. In our model, myogenic oscillations were phenomenologically considered as forcing tension applied to the vessel wall from VSMs. The model phenomenologically takes into account the modulation of SMC activity by vascular endothelium signaling (Koenigsberger et al., 2005; Delgado et al., 2010; Rapoport, 2014). One of the molecular mechanisms of this modulation may be explained through the release of vasoactive substances such as e.g., nitric oxide and/or Endotheli-1 expression as a result of activation of mechanosensitive channels in response to the shear stress (Vallance and Hingorani, 1999). Another mechanism of ECs and SMCs communication can be implemented here through calcium fluctuation in gap junction between SMCs and ECs. This mechanism was considered in the vasomotion model developed by Koenigsberger et al. (2005). The vasomotion model used in our paper [Equations (6) and (7)] phenomenologically describes ECs and SMCs communication through the negative feedback control in this coupled systems that is represented by negative terms in the right hand sides of these equations.

The coupled model of vasomotion which describes forced relaxation oscillation of active arterioles was solved numerically and spectral characteristics of its oscillatory solutions were calculated. The analysis of the power density spectrum of arteriole wall oscillations revealed two main frequency bands. The first one is in the frequency region of 0.01 Hz and corresponds to the relaxation oscillation in the coupled system. The second band in the range of 0.1 Hz is defined by forced myogenic oscillation. Analysis of the model showed that frequency patterns of vasomotion spectrum significantly depend on the amplitude of the myogenic oscillation. Amplitude of low-frequency limit-cycle oscillation in the range of 0.01 Hz decreases when amplitude of myogenic oscillations increases. This transformation of the spectrum at low- and highfrequency regions is defined by the phenomenon of subharmonic entrainment and locked mode in the forced relaxation oscillatory models, that was established well for coupled nonlinear models such as the forced van der Pol model (Jordan and Smith, 2007).

Note, in the model, we took into account only interaction between myogenic and endothelial oscillations and did not consider their coupling with neurogenic modes of vasomotions. The theoretical approach to consideration of this coupling in vasomotion models was discussed in Stefanovska (1999). Authors suggested that the coupling of the different modes ensures a coordinative control of the vascular function by different physiological mechanisms. Authors also mentioned that a study of these mode coupling can bear practical implications for diagnosis of microvascular circulation abnormalities (Stefanovska, 1999).

The transformation of the spectral density of vasomotion obtained in the model at the change of myogenic activity are in a satisfactory agreement with experimental data on the the LDF spectra measured in the stroke affected and unaffected hemispheres in AIS patients. A decrease in myogenic component in the frequency range of 0.1 Hz is accompanied by an increase

in endothelia component in the region of 0.01 Hz of the spectra. Base on this comparison, we assumed that declining in myogenic oscillation amplitude observed in the hemisphere affected by AIS is accompanied by an increasing amplitude of low frequency component in the frequency range of 0.01 Hz. Note that similar reciprocal relationship among different components of active oscillations in the LDF wavelet spectra was observed in LDF investigation of microcirculation (Krupatkin, 2009). Author assumed that disappearance of the myogenic components in the LDF wavelet spectrum is accompanied by increase in other components of the spectrum. The application of the vasomotion model (Arciero and Secomb, 2012) to analyse experimental data allowed us to suggest the possible mechanism of intraspectral changes of the LDF spectra measured in the patients with unilateral ischemic stroke.

In our theoretical description of the clinical data, we assumed identical parameters of the cerebrovascular systems of both hemispheres in healthy conditions. In the model, the asymmetry in the LDF signal arises from a relative difference in myogenic regulation of cerebrovascular dynamics in different hemispheres. Imbalance in endothelial signaling between hemispheres was not taken into account in the model. Although this factor related to endothelia dysfunction may affect asymmetry of the LDF signals in the endothelia frequency range (VI) of 0.005–0.0095 Hz. Another restriction of our theoretical approach to modeling of the LDF data is that we did not consider network vessel structure which forms the resultant LDF signal. In the simplified approach, we used the model of a single vessel with fixed pressure. To expand this model to the response of network vessel structure, multi-compartment approaches can be used which were developed to model vascular pathway (Ursino et al., 1996; Lanzarone et al., 2007; Arciero and Secomb, 2012).

A limitation of our experimental investigation is the absence of a control group for reliable estimation of inter-hemispheric variability (IHV) for AIS patients. This problem was also arisen in previous studies of the IHV and CA in LFO range of vasomotions (Phillip et al., 2012, 2013). In order to solve this problem, investigation of the IHV in healthy cohort as a control group was carried out using NIRS method (Phillip et al., 2012, 2013). Based on the absence of noticeable interhemispheric difference in the obtained data for the healthy cohort, authors suggested full synchronization of cerebrovascular regulation between hemispheres for healthy population and referred the observed IHV to impairment of perfusion and CA in affected and contralateral hemispheres in patients with occlusive disease (Phillip et al., 2013). Our data on the IHV in AIS patients conform in part with the IHV observation in LFO frequency region of vasomotion in the investigation of patients with AIS, carotid artery stenosis and cerebrovascular occlusion disease (COD) performed by NIRS and transcranial Doppler sonography (Reinhard et al., 2003; Phillip et al., 2013; Phillip and Schytz, 2014). In these studies, abrogation in inter-hemispheric synchronicity was defined mainly by the measurement of phase shift between ABP and oxyHb (Phillip et al., 2012; Phillip and Schytz, 2014) and between ABP and CBF (Reinhard et al., 2003) in the LFO range that showed impairment of CA at these cerebrovascular diseases. Amplitude ratio of ABP to oxyHb in the LFO range of spectrum was measured in COD (Phillip et al., 2012) and AIS patients (Phillip and Schytz, 2014) that showed higher ratio for the affected side compared to the contralateral one (ratio 1.3 ± 0.3 for AIS). While the measurement of the inter-hemispheric amplitude of oxyHb in the LFO range in AIS patients showed no significant difference in amplitude of the oxyHb LFOs on the strokeaffected side compared to the contralateral side (0.014 ± 0.03 vs. 0.017 ± 0.03; Phillip and Schytz, 2014). This result contradicts our data on declining spectral components in the frequency region near 0.1 Hz measured in stroke-affected hemisphere in 83.3% ischemic patients. Authors explained the absence of the amplitude difference in the stroke patients possible compensatory dilatation of arterioles distal to ischemic area as well as by LFOs decreasing with age and macroangiopathy in patient group (Phillip and Schytz, 2014). This discrepancy between our results and the results obtained in Phillip and Schytz (2014) requires further investigation of the inter-hemispheric variation with a wide range of patients by spectral characteristic monitoring in pre and in-treatment period up to complete or partial re-establishing of inter-hemispheric synchronization due to compensatory processes in affected and/or non-affected hemispheres.

# CONCLUSION

The preliminary results obtained in this pilot study showed that the inter-hemispheric variation in intra-spectral characteristics of LDF signals measured in AIS patients can be a prognostic factor for diagnostics of the cerebrovascular disease. Wavelet analysis of the LDF signals and theoretical investigation of vasomotions revealed that such prognostic factor may be the relationship between the spectrum components lying in the physiological relevant frequency bands related to myogenic and endothelial oscillations. Exploitation of the intra-spectral correlations in the inter-hemispheric variability of the LDF spectra can give useful information additionally to inter-spectral characteristics such as phase shift and correlation coefficients obtained in

### REFERENCES


NIRS and LDF investigation of cerebrovascular autoregulation impairments in patients with cerebrovascular diseases. For further development of this technique and translation of it to clinics, we suggest performing long term monitoring of the interhemispheric variation by LDF and NIRS and verify the obtained data by another group of patients with double blinded clinical trial. Spectral data collected in this monitoring will be valuable for estimation of population variants in the inter-hemispheric variability and reconciliation of the results obtained in different clinical trials. Analysis of these data along with computational modeling of vasomotion will help to understand physiological mechanism underlying inter-hemispheric variation in spectral characteristics at pre- and in-treatment period. The results of this analysis will facilitate design predictive and preventive intraand inter-spectral markers of the prevention, progression and treatment of cerebrovascular diseases.

# AUTHOR CONTRIBUTIONS

ER, VS, AK, and SS conceived main idea of the study; VS and AK designed the experiments; AA and MZ carried out a clinical investigation and data acquisition; AG conducted data analysis, modeling, and drafted the manuscript; SS, VS, and AK analyzed, interpreted the data, and contributed to the preparation of a final version of the manuscript.

#### ACKNOWLEDGMENTS

AG acknowledges Scottish Informatics and Computer Science Alliance (SICSA) for personal support and ER thanks for partial support Russian Ministry of Science and Education through grant No. 12.1223.2017/54

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2017.00160/full#supplementary-material

hemodynamics oscillations in TBI population. Brain Res. 1639, 194–199. doi: 10.1016/j.brainres.2016.02.018


patients with carotid stenosis. J. Cereb. Blood Flow Metab. 19, 460–465. doi: 10.1097/00004647-199904000-00012


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Goltsov, Anisimova, Zakharkina, Krupatkin, Sidorov, Sokolovski and Rafailov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Physical and Chemical Processes and the Morphofunctional Characteristics of Human Erythrocytes in Hyperglycaemia

Victor V. Revin<sup>1</sup> \*, Natalia A. Klenova<sup>2</sup> , Natalia V. Gromova<sup>1</sup> , Igor P. Grunyushkin<sup>1</sup> , Ilia N. Solomadin<sup>1</sup> , Alexander Y. Tychkov <sup>1</sup> , Anastasia A. Pestryakova<sup>2</sup> , Anna V. Sadykhova<sup>2</sup> , Elvira S. Revina<sup>1</sup> , Ksenia V. Prosnikova<sup>1</sup> , Jean-Christophe Bourdon<sup>3</sup> and Nikolai Zhelev <sup>4</sup> \*

*<sup>1</sup> Department of Biotechnology, Bioengineering and Biochemistry, Faculty of Biotechnology and Biology, Federal State-Financed Academic Institution of Higher Education, National Research Ogarev Mordovia State University, Saransk, Russia, <sup>2</sup> Department of Biochemistry, Biotechnology and Bioengineering, Faculty of Biology, Federal State-Funded Educational Institution of Higher Professional Education, Samara State University, Samara, Russia, <sup>3</sup> Division of Cancer Research, Ninewells Hospital and Medical School, University of Dundee, Dundee, United Kingdom, <sup>4</sup> CMCBR, Abertay University, Dundee, United Kingdom*

#### Edited by:

*Sergei G. Sokolovski, Aston University, United Kingdom*

#### Reviewed by:

*Natalia Pavlovna Akentieva, Institute of Problems of Chemical Physics (RAS), Russia Irma Borisovna Mosse, Institute of Genetics and Cytology at Belorussian Academy of Sciences, Belarus*

#### \*Correspondence:

*Victor V. Revin revinvv2010@yandex.ru Nikolai Zhelev n.zhelev@abertay.ac.uk*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *31 May 2017* Accepted: *07 August 2017* Published: *30 August 2017*

#### Citation:

*Revin VV, Klenova NA, Gromova NV, Grunyushkin IP, Solomadin IN, Tychkov AY, Pestryakova AA, Sadykhova AV, Revina ES, Prosnikova KV, Bourdon J-C and Zhelev N (2017) Physical and Chemical Processes and the Morphofunctional Characteristics of Human Erythrocytes in Hyperglycaemia. Front. Physiol. 8:606. doi: 10.3389/fphys.2017.00606* Background: This study examines the effect of graduated hyperglycaemia on the state and oxygen-binding ability of hemoglobin, the correlation of phospholipid fractions and their metabolites in the membrane, the activity of proteolytic enzymes and the morphofunctional state of erythrocytes.

Methods: Conformational changes in the molecule of hemoglobin were determined by Raman spectroscopy. The structure of the erythrocytes was analyzed using laser interference microscopy (LIM). To determine the activity of NADN-methemoglobinreductase, we used the P.G. Board method. The degree of glycosylation of the erythrocyte membranes was determined using a method previously described by Felkoren et al. Lipid extraction was performed using the Bligh and Dyer method. Detection of the phospholipids was performed using V. E. Vaskovsky method.

Results: Conditions of hyperglycaemia are characterized by a low affinity of hemoglobin to oxygen, which is manifested as a parallel decrease in the content of hemoglobin oxyform and the growth of deoxyform, methemoglobin and membrane-bound hemoglobin. The degree of glycosylation of membrane proteins and hemoglobin is high. For example, in the case of hyperglycaemia, erythrocytic membranes reduce the content of all phospholipid fractions with a simultaneous increase in lysoforms, free fatty acids and the diacylglycerol (DAG). Step wise hyperglycaemia in incubation medium and human erythrocytes results in an increased content of peptide components and general trypsin-like activity in the cytosol, with a simultaneous decreased activity of µ-calpain and caspase 3.

Conclusions: Metabolic disorders and damage of cell membranes during hyperglycaemia cause an increase in the population of echinocytes and spherocytes. The resulting disorders are accompanied with a high probability of intravascular haemolysis.

Keywords: human erythrocytes, hyperglycaemia, hemoglobin system, membrane lipids, proteolytic activity

**63**

# INTRODUCTION

Various physiological, pathological and nutritional conditions such as physical activity, large amounts of sweet food, emotional stress, metabolic syndrome, and diabetes are accompanied by high level of glucose in blood plasma. The high content of glucose in plasma accelerates the probability of non-enzymatic glycosylation of proteins, which induce damage to the cell membrane due to nonspecific aggregation of protein molecules and changes in protein-protein and protein-lipid interactions (Vasilyeva, 2005).

Taken together, these changes initiate the rapid aging of cells and the human organism. Metabolic syndrome significantly accelerates the development of atherosclerotic vascular damage and provokes earlier disability and death. During metabolic syndrome, which is currently the most common pathology of metabolic disorders, glycosylation of erythrocytic membrane proteins induces the impairment of rheological parameters of blood, low deformability and mobility of erythrocytes, high aggregation of erythrocytes and thrombocytes, high blood viscosity, and arterial hypertension (Shilov et al., 2008).

In addition, glycosylation of erythrocytic membrane proteins and hemoglobin during hyperglycaemia increases adhesion to endothelial cells, resulting in membrane destabilization (change in the asymmetry of membrane phospholipids), changes in viscoelastic properties of cells and their morphology (Riquelme et al., 2005). Taken together, these changes can impair the oxygentransport function of erythrocytes and reduce erythrocyte lifespan. Furthermore, the number of damaged circulating cells and, aging erythrocytes will increase (Lang et al., 2006; Mindukshev et al., 2010).

The biochemical mechanisms of impaired growth in human erythrocytes during the development of hyperglycaemia have not been sufficiently investigated. In particular, there are scarce data on the composition and status of the lipid phase of the membranes, the relationship of these processes with the activity of methemoglobin formation and the activity of apoptotic enzymes. Moreover, there is a lack of data in the literature on the effect of these processes on the morphofunctional state of erythrocytes and their oxygen-transport properties.

Therefore, we aimed to perform a comprehensive study of the effects of graduated hyperglycaemia on the composition of phospholipids, the activity of proteolytic enzymes, and, the consequent effect of ongoing processes on the morphofunctional state of erythrocytes and the oxygen-transport properties of hemoglobin.

# MATERIALS AND METHODS

#### Blood Samples

This study was performed using pure erythrocyte fractions isolated from freshly obtained, anticoagulated donor blood from the regional station of blood transfusion. The average age of the donors was 28.4 ± 1.7 years. The research was approved by the Local Ethics Board at Mordovia State University in accordance with the principles of Good Clinical Practice (protocol number 12 of 17 September 2014). Informed consent statements were signed by all donors participating in the experiment.

Erythrocytes were washed 4 times with PBS until pure fractions were formed and centrifuged at 600 g for 10 min at temperature +4 ◦C. Incubation of the erythrocytes was performed in a Ringer's solution at a ratio of 2:1 at 37◦C for 30 min with gentle mixing every 5 min (to avoid haemolysis). For incubation, a Ringer's solution was used containing as follow (in mM/l): NaCl–136.75; KCl–to 2.68; NaHCO3– 11.9 per; Na2HPO4–0.342; MgCl2–0.105; CaCl2–1.8, glucose–5 (normoglycaemia).

#### Hyperglycemia Model

Hyperglycaemic conditions were characterized by glucose content, in which the graduated concentration was accompanied by an equimolar decrease in the content of sodium chloride: 10 mM/l glucose–NaCl: 131.75 mM/l; 15 mM/l glucose–NaCl: 126.75 mm/l; 20 mM/l NaCl; 121.75 mM/L. The incubation medium was separated by centrifugation in the same manner, and erythrocytes were analyzed using spectrophotometry (Zavodnik and Lapshina, 1996) to determine the correlation of the hemoglobin forms.

### Conformation and Properties of Hemoglobin

Conformational changes in the hemoglobin molecule were determined by Raman scattering on In Via Renishaw (UK)



\**P*< *0.05 as related to normoglycemia;*

∆*P* < *0.05 as related to previous degree of hyperglycemia.*


microscope with a short-focus high-luminosity monochromator (focus distance–250 mm).

For excitation of the Raman spectra we used a laser with a wavelength of 532 nm, max radiation power of 100 mW, and 100 × objective. The data recorder was a –CCD detector (1,024 × 256 pixels with Peltier cooling to −70 ◦C) with a grating of 1800 lines/mm. The digitized spectra were processed in WIRE 3.3 (part of the unit's software). To analyse the conformation of globin haematoporphyrin, we used specific bands of the RS spectrum, which enabled an estimate of the relative amount of oxyhemoglobin and the ability of hemoglobin to bind to and isolate ligands, including oxygen (Brazhe et al., 2009).

#### Laser Interference Microscopy

The structure of the erythrocytes was analyzed using laser interference microscopy (LIM) in vitro (Byazhe et al., 2006; Yusipovich et al., 2011; Revin et al., 2016) using MII-4 sysytem (Russia). The measurements were performed at room temperature, and a suspension of erythrocytes in the incubation medium (1:2) was placed on mirror glass. The smear was prepared and covered with cover glass. Images of 10 sites, with a monolayer arrangement of cells in an interference channel, were obtained, with reflected light in each sample. The images were processed using FIJI (Schindelin et al., 2012). The structure of the erythrocytes was assessed by registering the average value of the optical path difference (OPD) and phase image area using at least 100 cells from each sample. The phase volume of the erythrocyte was calculated using the following formula:

$$\mathbf{V}\_{\text{cell}} = \mathbf{F}\_{\text{mean}} \cdot \mathbf{S} / \mathbf{n}\_{\text{cell}} - \mathbf{n}\_{\text{m}}$$

Where Fmean is the mean value of the optical path difference, proportional to the thickness of erythrocyte; S is the phase image area of the cells; ncell is the refractive index of erythrocyte, equal to 1.405; n <sup>m</sup> is the refractive index of the surrounding solution (1.333).

#### Determination of the Activity of NADN-Methemoglobinreductase

To determine the activity of NADN-methemoglobinreductase, we used the P.G. Board method (Board, 1981). The degree of glycosylation of the erythrocyte membranes was determined using a method previously described by Felkoren et al. (1991) . Before defining the proteolytic enzymes, the suspension of erythrocytes was haemolysed by the addition of a buffer of 20 mM three-HCl containing 2 mM EDTA, pH 7.5 in a ratio of 1:9. The haemolysate was maintained at a temperature 2-4 ◦C for 15 min and centrifuged at 16,000 × g for 40 min.

The supernatant was used to determine the activity of µ-calpain with the release of enzyme using ion exchange chromatography (column 3 × 15; DEAE-cellulose) and eluted by a gradient of 0.1–0.4 M NaCl. Next, the mean calpain activity of the fractions was determined using incubation medium previously described by Sorimachi et al. (1997) (imidasole buffer, 4% casein, 50 mM of CaCl <sup>2</sup>, 50 mM of cysteine) (Stroev et al., 1991; Elce John, 2000; Sorimachi et al., 2000). Then, the general proteolytic activity, in ascending order of the absorption level at a

TABLE

2


 hemoglobin

 in

hyperglycaemia.



\**P*< *0.01 as related to normoglycemia;*

◦*P*< *0.01 as related to previous level.*



\**P* < *0.01 as related to normoglycemia;*

◦*P* < *0.01 as related to the level of 10 mM/l.*

wavelength of 280 nm, was measured following incubation of the haemolysate and protein deposition by 5% trichloroacetic acid (TCA) (Bazarnova et al., 2008).

The release of fractions with calpain activity occurred during movement in a reverse gradient from 0.2 to 0.1 M of sodium chloride immediately after the release of hemoglobin. The activity of µ-calpain was calculated as the difference between activity with and without the inhibitor in the incubation medium (phenylmetylsulfonylfluoride, 2 mM PMSF, Sigma, USA). The content of peptides in the incubation medium and erythrocytes was determined using the Lowry method, with the Bio-Rad DC Protein Assay a set of protein assay Dc reagents (Bio-Rad). The content of peptides in the erythrocytes was determined after deproteinisation by 5% of TCA. The active concentration of caspase-3 in erythrocytes was recorded using an enzyme immunoassay (BD Biosciences, USA) with Stat Fax 3200 microplate reader (USA).

#### Analysis of Membrane Phospholipids

To analyse the state of the membrane phospholipids, we isolated membranes from the haemolysate using 5 mM NaH2PO4+ 0.5 mM PMSF (phenylmetylsulfonylfluoride) solution, cooled to 0 ◦C, pH 8.0 in a ratio of 1:20. The mixture was incubated for 10 min at 4◦C and then centrifuged at 20,000 × g for 40 min (0◦C). The supernatant was removed and the residue was resuspended in lysis solution and centrifuged in the same manner. The sample was washed three times. Lipid extraction was performed using the Bligh and Dyer method (Bligh and Dyer, 1959).

To separate the phospholipid fractions we used onedimensional chromatography combined with chloroform/ TABLE 5 | Content of protein and peptide compounds in the incubation medium and peptides in erythrocytes (during normo- and hyperglycaemia).


\**P*< *0.01 as related to normoglycemia;*

◦*P*< *0.05 as related to normoglycemia.*

methanol/glacial acetic acid/water in a ratio of 60/50/1/4 (Evans et al., 1990). Chromatographic separation was performed in a thin layer of silica gel deposited on a glass plate. Standard plates from HPTLC Silicagel 60 F254 (Merck. Germany) were used. To separate DAG and FFA we used the combination of heptane/diethyl ether/glacial acetic acid (60/40/2 by volume). Detection of the phospholipids was performed using V. E. Vaskovsky method (Vaskovsky et al., 1975).

The amount of phospholipid fractions and free fatty acids in the erythrocytic membranes was determined using a densitometer TCL Scanner 3 (Gamag, Switzerland) at an absorption wavelength of 360 nm using a deuterium lamp and winCATC software.

#### Statistical Analysis

The statistical processing was performed using the Statistika 0.06 software package. To determine the significance levels, the Student's test was used. Repeats in the variation series of different indicators were from 8 to 20 values.

#### RESULTS

Hyperglycaemic conditions are characterized by a low affinity of hemoglobin to oxygen, which is manifested as a parallel decrease in the content of hemoglobin oxyform and the growth of desoxyform (**Table 1**).

#### TABLE 6 | Activity of erythrocyte enzymes during graduated hyperglycaemia.


\**P* < *0.01 as related to normoglycemia;*

◦*P* < *0.05 as related to normoglycemia.*

TABLE 7 | Morphological parameters of erythrocytes in normo- and hyperglycaemia, according to LIM results.


\**P*< *0.01 as related to normoglycemia;*

∆*P*< *0.01 as related to previous level of hyperglycemia.*

Upon a 30-min incubation of the cells, an increase in the glycosylation of hemoglobin occurred only when the glucose concentration was 20 mM/l, whereas a small but reliable growth in methemoglobin formation, increase in membrane-bound form of hemoglobin and degree of glycosylation of proteins in erythrocytic membranes has been previously observed when the glucose level exceeded a two-fold level (10 mm/l) (**Table 1**).

The data obtained using the spectrometric measurements were confirmed by Raman scattering (**Table 2**).

We found a reliable decrease in the relative content of oxygenated hemoglobin, with an average of 2%. Moreover, we obtained similar spectrophotometric percentages of hemoglobin forms (**Tables 1**, **2**). The conformational changes observed indicate an increased ability of hemoglobin to bind to ligands and the reliable occurrence of low hemoglobin affinity to oxygen in the case of 20-mM/l hyperglycaemia, at approximately 11%. Furthermore, there is a low intensity of symmetric and asymmetric vibrations of pyrrole rings, which indicated a less conformational mobility of the haeme structures and an impairment of the effective binding of ligands (**Table 2**).

At 15- and 20-mM/l glucose hyperglycaemia, we detected a high level of lysophosphatidylcholine (LFH) in the erythrocyte membrane, and the content of phosphatidylcholine (PFS) was reliably low under in conditions of increasing amounts of glucose, i.e., up to 10 mM/l. In addition, we recorded a low level of sphingomyelin (SM) and fractions of phosphatidylinositol + phosphatidylserine (PI + PS). Strong hyperglycaemia (20 mm/l) is characterized by a low content of all fractions with a simultaneous increase of in LPC (**Table 3**).

The percentage of increased lysophosphatides under strong 20-mm/l hyperglycaemia was 160%. The SM level was 2 to 3 times lower, with PC–3 times; PI+PS–2.35%; and PEA–4 times. The low content of phospholipids and high content of lysophosphatides indicated the activation of phospholipase A and C, which was accompanied by loosening of the membrane and higher ion permeability (Mills and Needham, 2005).

Analysis of the content of free fatty acids (FFA) and the product of the reaction catalyzed by phospholipase C and diacylglycerol (DAG) showed growth that was proportional to the level of hyperglycaemia, and the percentage of growth with 20 mM/l hyperglycaemia is 133 and 182%, respectively (**Table 4**).

These results clearly demonstrate the activity of phospholipase.

The normal correlation of phospholipid fractions predetermines the effective regulation of active and passive transport of substances, the sensitivity of cells to the action of ligands, and the activity of membrane-bound enzymatic systems. The optimal operation of ion-carrying systems (Ca2<sup>+</sup> and Na+, K+- ATP-PS) is facilitated by stable PEA content. During severe hyperglycaemia (20 mm/l), a sharp decrease in PEA content will be accompanied by a defect in ion-carrying processes (**Table 3**). Moreover, it will trigger a general change in the structure of the membrane phospholipid bilayer, as PEA is a structural phospholipid (Delaunay, 2002).

Importantly, the low content of PEA causes a disorder in endoglobular homeostasis and inhibition of the antioxidant activity of cells in (Afanasiev et al., 2007). The low content of phospholipids with polyunsaturated fatty acids (PC, PI, PEA) triggers defatting of membranes and a high correlation of cholesterol/PL, which is accompanied by changes in physicochemical properties, namely, high microviscosity. These findings indicate that hyperglycaemia is followed by a disruption in membrane permeability.

In hyperglycaemia, human erythrocytes demonstrate a high content of peptide components, which was specifically high when the content of glucose in the incubation medium was 10 mm/l (**Table 5**). A further increase in the level (degree) of hyperglycaemia was accompanied by an increased release of protein compounds into the incubation medium due to higher cell membrane permeability (**Table 5**).

Measurement of the proteolytic activity revealed interesting patterns. The release of µ-calpain fractions and identification of the average activity of the enzyme showed a decrease with a parallel increase in the degree of hyperglycaemia, which was statistically significant when the concentration of glucose reached 20 mM/l. A low activity of caspase 3 was detected, which was also strongest under conditions of severe hyperglycaemia (**Table 6**). Increased proteolytic activity can be detected only by the identification of the overall trypsin-like activity of cytosol with respect to the high content of peptides in the haemolysate, without the isolation of enzymes and addition of third proteolytic substrates (**Table 6**).

Laser interference microscopy showed distinct structural changes of erythrocytes with increased concentrations of glucose in the incubation medium (**Table 7**, **Figures 1**–**4**).

During normoglycaemia, erythrocytes appeared similar to discocytes, with a normal distribution of hemoglobin inside the cells. During hyperglycaemia, the estimated parameters were high overall. In addition to the relative areas of the erythrocytes, all indicators were high in conditions of 20-mm/l hyperglycaemia (**Table 7**).

With a high content of glucose in the incubation medium, echinocytes accumulated in the population of cells, and the cell profile demonstrated visible outgrowths on the erythrocyte surface. At up to 15-mm/l glucose, spherocytes appeared, and the cells appeared more swollen (**Figure 3**). At 20-mm/l glucose, the cells were slightly compressed compared with 15-mm/l glucose; however, other changes became more robust (**Table 7**, **Figure 4**).

### DISCUSSION

Disorder in the metabolic and functional state of human erythrocyte membranes at a high concentration of glucose in the incubation medium occurs primarily due to the general activation of biochemical and physicochemical processes in cells. The shift in the dissociation curve in the direction of deoxygenation causes increased methemoglobin formation as DOHb, which is less resistant to auto-oxidation compared tooxyform (Zavodnik and Lapshina, 1996; Ivanov, 2001). The gradual reduction of the conformational mobility of haeme structures (**Table 2**) is most likely caused by the high degree of the glycosylated protein component of hemoglobin, its permeability in the membrane, followed by a gradual increase in membranebound hemoglobin (**Table 1**).

The observed increase in the degree of glycosylation of membrane proteins induces changes in the state of membrane phospholipids. The change in protein-lipid interactions results in the activation of membrane phospholipase, causing an increase in LPC and DAG, and results in a high content of FFA in erythrocytic membranes. Membrane permeability under these conditions increases, and the peptide and protein compounds are released into the incubation medium. Some reduction in their release at 20-mm/l hyperglycaemia may be accounted for by the low activity of µ-calpain and caspase 3 due to less affinity to the glycosylated form of hemoglobin. However, their content in the erythrocytes and incubation medium was higher compared to normoglycaemia, most likely due to the continued high trypsinlike activity of the cell cytosol.

Metabolic disorders and injury to the cell membrane trigger morphological changes in erythrocytes. Previous research findings have shown a gradual accumulation of structurally damaged cells in the population. First, in 10-mM/l hyperglycaemia, the number of echinocytes increases (reversible morphological changes) and the phase volume and thickness of the cells increases. Next, stomacytes and spherocytes appear in the population. LIM microscopic results indicated strong cell swelling. In 20-mm/l hyperglycaemia, the number of spherocytes increased. Furthermore, the cells showed some shrinkage, but the ratio of surface/volume values reached a critical level. Thus, haemoglytic destruction is very likely.

The destruction of erythrocytes in the bloodstream in hyperglycaemia will result in a chronic inflammation response, which can damage the body's blood vessels.

Changes in the erythrocytic membrane not only result in morphological disorders of the whole erythrocyte, but they can also consequently result in conformational changes in hemoglobin as well as its oxygen-binding and oxygen-transport function. The low content of oxyhemoglobin, weak symmetric and asymmetric oscillations of pyrrole rings and increase in hemoglobin affinity to ligands, including oxygen, has been demonstrated.

# CONCLUSION

Taken together, we demonstrate the significant risk of the prolonged availability of high glucose concentrations in plasma. Nonspecific glycosylation of membrane proteins and erythrocytesin hemoglobin results in a weak affinity of hemoglobin to oxygen and its loss by cells enroute to tissues in the human body. Moreover, the resulting damage to membranes and cell metabolism increase the probability of an accumulation of functionally defective aging erythrocytes in the circulating population. More rapid spherization of erythrocytes, in the absence of physiological mechanisms of apoptosis, introduces the possibility of necrotic cell death in the bloodstream, resulting in the gradual development of chronic states of hypoxia and inflammation.

#### AUTHOR CONTRIBUTIONS

VR: designed experiments, interpreted data, wrote manuscript. NK, NG, IG, IS, AT, AP, AS, and ER: designed and performed experiments wrote manuscript. KP: performed experiments and

#### REFERENCES


analyzed results. NZ and JCB: made substantial contribution to the analysis and interpretation of the data, revised and critically reviewed the manuscript.

#### FUNDING

The authors are grateful for the support of the Russian Science Foundation, Grant no. 15-15-10025.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Revin, Klenova, Gromova, Grunyushkin, Solomadin, Tychkov, Pestryakova, Sadykhova, Revina, Prosnikova, Bourdon and Zhelev. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Optoacoustic Monitoring of Physiologic Variables

#### Rinat O. Esenaliev\*

*Laboratory for Optical Sensing and Monitoring, Department of Neuroscience and Cell Biology, Department of Anesthesiology, Center for Biomedical Engineering, University of Texas Medical Branch, Galveston, TX, United States*

Optoacoustic (photoacoustic) technique is a novel diagnostic platform that can be used for noninvasive measurements of physiologic variables, functional imaging, and hemodynamic monitoring. This technique is based on generation and time-resolved detection of optoacoustic (thermoelastic) waves generated in tissue by short optical pulses. This provides probing of tissues and individual blood vessels with high optical contrast and ultrasound spatial resolution. Because the optoacoustic waves carry information on tissue optical and thermophysical properties, detection, and analysis of the optoacoustic waves allow for measurements of physiologic variables with high accuracy and specificity. We proposed to use the optoacoustic technique for monitoring of a number of important physiologic variables including temperature, thermal coagulation, freezing, concentration of molecular dyes, nanoparticles, oxygenation, and hemoglobin concentration. In this review we present origin of contrast and high spatial resolution in these measurements performed with optoacoustic systems developed and built by our group. We summarize data obtained *in vitro*, in experimental animals, and in humans on monitoring of these physiologic variables. Our data indicate that the optoacoustic technology may be used for monitoring of cerebral blood oxygenation in patients with traumatic brain injury and in neonatal patients, central venous oxygenation monitoring, total hemoglobin concentration monitoring, hematoma detection and characterization, monitoring of temperature, and coagulation and freezing boundaries during thermotherapy.

#### Edited by:

*Alexey Goltsov, Abertay University, United Kingdom*

#### Reviewed by:

*Valeriy G. Andreev, Moscow State University, Russia Alexey Popov, University of Oulu, Finland*

#### \*Correspondence:

*Rinat O. Esenaliev riesenal@utmb.edu*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *19 August 2017* Accepted: *28 November 2017* Published: *12 December 2017*

#### Citation:

*Esenaliev RO (2017) Optoacoustic Monitoring of Physiologic Variables. Front. Physiol. 8:1030. doi: 10.3389/fphys.2017.01030* Keywords: optoacoustic, photoacoustic, monitoring, imaging, sensing, physiologic

# INTRODUCTION

Optical techniques can be used for noninvasive measurements of physiologic variables, functional imaging, and hemodynamics monitoring. These measurements are mostly based on high optical absorption contrast of tissue chromophores (Welch and van Gemert, 2011; Tuchin, 2016). However, optical techniques have drawbacks associated with limited resolution due to strong light scattering in tissues.

Optoacoustic technique utilizes thermoelastic generation of optoacoustic (ultrasound) waves in tissues by short optical pulses. Optoacoustic pressure wave amplitude is linearly dependent on the absorption coefficient. Time-resolved detection of the optoacoustic waves yields high (optical) contrast and high (ultrasound) resolution that can be used for imaging, sensing, and monitoring in tissues.

**71**

Since early 1990s we have been working on biomedical optoacoustics and proposed to use it for many diagnostic applications, developed and built optoacoustic systems, and tested them in tissues phantoms, tissues in vitro and in vivo, animal models, and clinical studies (Esenaliev et al., 1993, 1996, 1997, 1998a, 1999b, 2002a, 2004a,b; Oraevsky et al., 1998; Larin et al., 2001, 2002, 2005; Larina et al., 2005; Petrova et al., 2005, 2009; Petrov et al., 2005, 2006, 2012a,b, 2014, 2016, 2017a,b; Brecht et al., 2007; Patrikeev et al., 2007; Herrmann et al., 2017). In this review we present origin of contrast and high spatial resolution in physiologic measurements and discuss results obtained with the optoacoustic systems in vitro, in animals, and in humans on monitoring of physiologic variables.

## ORIGIN OF HIGH CONTRAST AND RESOLUTION IN OPTOACOUSTIC MEASUREMENTS

#### Theoretical Background

The thermoelastic mechanism of optoacoustic wave generation is based on absorption of light energy in a medium followed by temperature rise and thermal expansion in the medium. The thermal expansion of the irradiated medium induces mechanical stress (pressure rise). A short optical pulse with the incident fluence, Fo, induces a pressure rise, P(z), in the medium upon condition of stress confinement:

$$P(z) = (\beta c\_s^{-2}/C\_p)\mu\_a F = \Gamma \mu\_a F(z) = \Gamma \mu\_a F\_o \exp(-\mu\_a z) \tag{1}$$

where β [1/◦C] is the thermal expansion coefficient; c<sup>s</sup> [cm/s] is the speed of sound; C<sup>p</sup> [J/g◦C] is the heat capacity at constant pressure; F(z) [J/cm<sup>2</sup> ] is the fluence of the optical pulse; and µ<sup>a</sup> [cm−<sup>1</sup> ] is the absorption coefficient of the medium. The generated optoacoustic pressure can be expressed in J/cm<sup>3</sup> or in bar (1 J/cm<sup>3</sup> = 10 bar). The combination of the thermophysical parameters, βc 2 s /C<sup>p</sup> in Equation (1) represents the Grüneisen parameter, Ŵ (dimensionless). The exponential light attenuation in the medium is represented by exp(−µaz).

The Equation (1) is valid upon the stress-confinement condition when pressure relaxation is negligible during the heat deposition. The stress-confined condition is satisfied when light pulse duration, τp, is shorter than the stress relaxation time in the irradiated volume,τstr:

$$
\pi\_p < \pi\_{str} = \frac{1}{\mu\_a c\_s} \tag{2}
$$

Nanosecond pulses can be used to generate conditions of stress confinement for many biomedical optoacoustic applications including monitoring of physiologic variables (Esenaliev et al., 1993, 1996, 1998a; Larin et al., 2001, 2002, 2005; Larina et al., 2005; Petrova et al., 2005, 2009; Petrov et al., 2005, 2006; Brecht et al., 2007).

According to Equation (1), optoacoustic pressure amplitude is proportional to the Grüneisen parameter, fluence, and absorption coefficient of the medium, while the pressure spatial profile is dependent on the absorption coefficient. Since z and t are related by the simple equation:

$$z = c\_t t \tag{3}$$

the spatial distribution of optoacoustic pressure P(z) is detected by a wide-band acoustic transducer as a temporal profile P(t):

$$P(t) = \Gamma \mu\_a F\_o \exp(-\mu\_a c\_s t) \tag{4}$$

Therefore, by analyzing the amplitude and/or temporal profile of optoacoustic waves, one can measure the absolute value of the absorption coefficient of the irradiated medium.

Most tissues are strongly scattering media in the optical spectral range. In addition to the absorption coefficient, two other major optical parameters are responsible for distribution of light in tissues: scattering (µs) and effective attenuation (µeff) coefficients. Attenuation of diffusively scattered light depends on the effective attenuation coefficient which is related to µa, µ<sup>s</sup> , and the anisotropy factor (g) as:

$$
\mu\_{\rm eff} = \{ \Im \mu\_a [\mu\_a + \mu\_s (1 - \mathfrak{g})] \}^{1/2} \tag{5}
$$

where µs(1–g) is the reduced scattering coefficient, µ<sup>s</sup> ′ (Welch and van Gemert, 2011). Light penetration depth in tissues is defined as 1/µeff. Distribution of laser fluence and, therefore, pressure in tissue (not very close to the surface) is dependent on optical absorption and effective attenuation coefficients:

$$P(z) = \Gamma \mu\_d k F\_o e^{(-\mu\_{\epsilon \mathcal{Y}} z)} \tag{6}$$

where k is the parameter resulted from multiple scattering in tissue and depends on absorption and scattering coefficients (Welch and van Gemert, 2011).

Absorption and reduced scattering coefficients of tissues are low in the near-IR spectral range (from 700 to 1,300 nm), that results in deeper penetration of near-IR radiation compared with that of other parts of the spectrum. Application of near-IR radiation allows for sufficient penetration of light in tissues for optoacoustic measurements of physiologic variables.

#### Monitoring of Temperature

Accurate temperature mapping with sub-mm spatial resolution may provide precise thermotherapy of abnormal tissues with minimal damage to surrounding normal tissues. Amplitude of optoacoustic pressure waves induced in many tissues increases with temperature, mostly due to temperature dependence of the thermal expansion coefficient (Esenaliev et al., 1998b; Larin et al., 2002, 2005; Larina et al., 2005). The temperature-dependent Grüneisen parameter can be expressed with an equation:

$$
\Gamma = A + BT \tag{7}
$$

where A and B are constants and T is temperature. One can modify (Equation 1) for the case of absorbing media without scattering as:

$$P(z) = \left(A + BT(z)\right)\mu\_a F\_o e^{(-\mu\_a z)}\tag{8}$$

and for the case of strongly scattering media in deeper (not in the subsurface) areas as:

$$P(z) = \left(A + BT(z)\right)k\mu\_a F\_o e^{(-\mu\_{\epsilon\mathcal{J}}z)}\tag{9}$$

where T(z) is the temperature distribution in tissue. One can rearrange (Equation 9) to obtain temperature distribution in tissue:

$$T(z) = C + D \cdot P(z) / P(z)\_{\Gamma = \text{To}} \tag{10}$$

where P(z)T=To is the optoacoustic pressure profile recorded at the initial temperature T<sup>o</sup> and C and D are parameters that are dependent on tissue properties. Therefore, by recording and analyzing the temporal optoacoustic pressure profile, one can reconstruct distribution of temperature during hyperthermia.

We experimentally demonstrated linear increase with temperature of optoacoustic pressure amplitude in tissue phantoms and tissues such as liver and myocardium (Esenaliev et al., 1998b). In another set of experiments, using rapid heating, we produced temperature gradients in tissue and tissue-like sample (Larina et al., 2005). During the heating we monitored the temperature distribution with optoacoustic system, while a multi-sensor temperature probe inserted in the samples measured actual temperature distribution. These studies demonstrated that the accuracy of optoacoustic temperature was better than 1◦C at the spatial resolution less than 1 mm.

#### Monitoring of Coagulation and Freezing Front

Tissue optical and thermophysical properties change due to coagulation or freezing. This may provide fast and accurate optoacoustic feedback during thermotherapy with heating or cooling agents. Because the optoacoustic wave amplitude and temporal parameters are dependent on tissue properties (Equation 6), detection and analysis of the optoacoustic waves during thermotherapy may be used for real-time monitoring of the extent of tissue coagulation or freezing with high resolution and contrast (Esenaliev et al., 1998a; Larin et al., 2002, 2005; Larina et al., 2005).

We performed high-resolution, real-time optoacoustic monitoring of tissue coagulation during conductive heating (Larina et al., 2005) or interstitial heating by CW laser light (Larin et al., 2005). Analysis of optoacoustic signal slopes was used for monitoring tissue heating and dimensions of coagulation lesions. Coagulation was induced in liver, myocardium, and prostate by interstitial CW Nd:YAG laser irradiation of the samples or by conductive heating. The optical properties did not change up to the coagulation temperature (about 53◦C), but sharply increased during heating up to 70◦C. The interstitial coagulation front was monitored in freshly excised canine tissues in real time with spatial resolution of about 0.6 mm. These results suggested that this technique may be used for real-time precise thermotherapy of malignant and benign lesions at depths of the order of centimeter.

Real-time monitoring of cooling and freezing of tissues, cells, and other biological objects with a high spatial and temporal resolution is necessary for selective cryoablation of cancer and benign tumors and for organs, tissues, and other biological objects in medicine and cryobiology. Using liquid nitrogen, we demonstrated that tissue hypothermia and freezing can be monitored with the optoacoustic technique because both amplitude and profile of the optoacoustic waves change with temperature during cooling and freezing (Larin et al., 2002). Sharp increase of the optoacoustic signal slope was detected between −2 and −4 ◦C resulted from the formation of the frozen zone in liver tissue. High spatial resolution (better than 0.5 mm) was obtained in these studies. Such resolution is sufficient for cryoablation monitoring with high precision.

TABLE 1 | Physiologic variables that were monitored in our optoacoustics works, origin of contrast, range, accuracy, precision of monitoring, and corresponding references.


# Monitoring of Exogenous Dyes and Nanoparticles

Optical absorption contrast can be used for monitoring in tissues of exogenous dyes and nanoparticles with high temporal and spatial resolution. Indocyanine green (ICG), an FDAapproved dye for intravenous injections has a high absorption in the near IR spectral range with a maximum at 800–805 nm. Optoacoustic measurements of ICG-produced signal can be useful for monitoring of cardiac output (CO), cardiac index (CI), blood volume (BV), and hepatic function. We measured amplitude and peak-to-peak amplitude of optoacoustic signals induced in whole arterial blood in vitro at clinically relevant ICG concentrations. Both amplitude and peak-to-peak amplitude of the signals linearly increased with ICG concentration with high correlation coefficients (R <sup>2</sup> = 0.990 and 0.991, respectively) (Prough et al., 2008). The blood effective attenuation coefficient derived from the optoacoustic signal slopes also increased linearly with ICG concentration with high correlation coefficient of R <sup>2</sup> = 0.986.

Nano- and microparticles that strongly absorb light can be used for photothermal therapy of abnormal tissues including tumors (Esenaliev, 1999a, 2000, 2016a). High sensitivity of the optoacoustic technique to absorption changes provides basis for monitoring of nanoparticle delivery into tumors and for tumor coagulation. We used the optoacoustic technique for real-time monitoring of nanoparticle accumulation in human tumors of nude mice and the nanoparticle-induced laser thermotherapy of these tumors (Esenaliev et al., 2007).

#### Monitoring of Hemoglobin Concentration

Hemoglobin has a high absorption coefficient in the visible and near-IR spectral range (Welch and van Gemert, 2011). We experimentally demonstrated that both the amplitude and spatial distribution of the generated optoacoustic pressure induced in blood are dependent on total hemoglobin concentration [THb] (Esenaliev et al., 2004a,b; Petrova et al., 2005).

High z-axial resolution of the optoacoustic technique permits direct [THb] measurements in blood vessels because the optoacoustic waves induced in blood arrive at the acoustic transducer at the time defined by Equation (3). Since the hemoglobin absorption coefficient is dependent on hemoglobin saturation (i.e., oxygenation), we used the wavelength of 805 nm (isobestic point) where oxy- and deoxygenated hemoglobin have same absorption (Esenaliev et al., 2004a,b; Petrova et al., 2005). This allows for [THb] measurement at any oxygenation.

### Monitoring of Oxygenation

One of the most important optoacoustic applications is oxygenation imaging, monitoring, and sensing (Esenaliev et al., 2002a,b). They can be used for diagnostics and management of large populations of patients including those with traumatic brain injury (TBI), circulatory shock, stroke, patients undergoing surgery, anemic patients, neonatal patients, and fetuses during late-stage labor (Petrova et al., 2005, 2009; Petrov et al., 2005, 2006, 2012a,b, 2014, 2016, 2017a,b; Brecht et al., 2007; Esenaliev et al., 2016b; Esenaliev, 2017; Herrmann et al., 2017).

We proposed and developed a noninvasive, optoacoustic diagnostic platform for measurements of oxygenation, hemoglobin concentration, and other important physiological parameters in tissues and specific blood vessels (Petrova et al., 2005, 2009; Petrov et al., 2005, 2006, 2012a,b, 2014, 2016, 2017a,b; Brecht et al., 2007; Prough et al., 2008; Herrmann et al., 2017). Because hemoglobin is a major chromophore in the near IR spectral range and its absorption depends on oxygenation, optoacoustics is suitable for monitoring of these physiologic variables.

# DISCUSSION

In early 1990s we performed first experimental studies on biomedical optoacoustics in tissues and mid 1990s obtained first results on optoacoustic microscopy and demonstrated optoacoustic signal detection in tissues at several centimeters depths that is well beyond the optical diffusion limit (five times greater than the light penetration depth defined as 1/µeff). High-resolution optoacoustic images of tissue phantoms were reconstructed by our group in early 2000s; since then biomedical optoacoustics/photoacoustics has grown tremendously. Recent original papers and reviews by other groups (and references therein) demonstrate remarkable progress in optoacoustic imaging, cancer detection, microscopy, functional imaging in small animals, optoacoustic instrumentation, dual modality optoacoustic/ultrasound imaging, contrast agents development, and other applications (Jaeger et al., 2013; O'Donnell et al., 2013; Su et al., 2013; Daoudi et al., 2014; Ellwood et al., 2014; Yao and Wang, 2014; Taruttis et al., 2015; Bourantas et al., 2016; Cai et al., 2016; Choi et al., 2016).

The results of our studies demonstrated high contrast and high resolution in optoacoustic monitoring, imaging, and sensing of the physiologic variables. The high contrast originates from thermophysical properties (Grüneisen parameter) and/or from optical properties (absorption and effective attenuation coefficient) of tissue. **Table 1** shows the physiologic variables that were monitored in our optoacoustic works, origin of contrast, range, accuracy, and precision of monitoring. Combination of the high contrast with high resolution (from microns to sub-mm) yields a promising diagnostic modality that can be used in clinics and biomedical research. Although light attenuation decreases the monitoring accuracy in thick tissues, high signal-to-noise ratio of optoacoustic signals allowed for accurate measurements of these physiologic variables at depth greater than the optical diffusion limit (**Table 1** and references therein).

# CONCLUSION

We proposed to use optoacoustics for imaging, monitoring, and measurements of a number of important physiologic variables and developed optoacoustic systems for these applications. The systems were successfully tested in tissue phantoms, tissues, animals, and human subjects. The obtained data suggest that this technology may be applicable to large populations of patients. We plan to further develop these systems for single measurement, continuous measurement, and monitoring, as well as for 2D, 3D, and 4D imaging of the physiologic variables.

## AUTHOR CONTRIBUTIONS

The author confirms being the sole contributor of this work and approved it for publication.

## ACKNOWLEDGMENTS

The author thanks Drs. Donald S. Prough, Yuriy Petrov, Irene Y. Petrov, Claudia S. Robertson, C. Joan Richardson, other members of the Department of Anesthesiology, Department of

# REFERENCES


Pediatrics and Neonatal Intensive Care Unit at UTMB, and the Department of Neurosurgery and Neurointensive Care Unit of Baylor College of Medicine. Grant support: NIH grants #R01EB00763 and #U54EB007954 from the National Institute of Biomedical Imaging and Bioengineering, #R01NS044345 and #R21NS40531 from the National Institute of Neurological Disorders and Stroke, #R41HD076568 and #R43HD075551 form the Eunice Kennedy Shriver National Institute of Child Health & Human Development, and #R41HL10309501 from the National Heart, Lung and Blood Institute, contracts from Noninvasix, Inc., the Texas Emerging Technology Fund, and the Moody Center for Brain and Spinal Cord Injury Research/Mission Connect of UTMB. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIBIB or NIH.


**Conflict of Interest Statement:** RE is a co-owner of Noninvasix, Inc., a UTMBbased startup that has licensed the rights to optoacoustic technology.

Copyright © 2017 Esenaliev. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Use of Microwave Radiometry to Monitor Thermal Denaturation of Albumin

Yuri Ivanov <sup>1</sup> \*, Andrey F. Kozlov <sup>1</sup> , Rafael A. Galiullin<sup>1</sup> , Vadim Y. Tatur <sup>2</sup> , Vadim S. Ziborov <sup>3</sup> , Nina D. Ivanova<sup>4</sup> , Tatyana O. Pleshakova<sup>1</sup> , Sergey G. Vesnin5,6 and Igor Goryanin7,8,9 \*

*1 Institute of Biomedical Chemistry, Moscow, Russia, <sup>2</sup> Foundation of Advanced Technologies and Innovations, Moscow, Russia, <sup>3</sup> Joint Institute for High Temperatures of Russian Academy of Sciences (RAS), Moscow, Russia, <sup>4</sup> Moscow State Academy of Veterinary Medicine and Biotechnology, Moscow, Russia, <sup>5</sup> RES LTD, Moscow, Russia, <sup>6</sup> Medical MicroWave Radiometry (MMWR) LTD, Edinburgh, United Kingdom, <sup>7</sup> School of Informatics, University of Edinburgh, Edinburgh, United Kingdom, <sup>8</sup> Biological Systems Unit, Okinawa Institute of Science and Technology, Okinawa, Japan, <sup>9</sup> Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences, Tianjin, China*

#### Edited by:

*Alexey Goltsov, Abertay University, United Kingdom*

#### Reviewed by:

*Julian Borejdo, Unversity of North Texas Health Science Center, United States Vladimir Danichev, Kurchatov Institute, Russia*

#### \*Correspondence:

*Yuri Ivanov yurii.ivanov@rambler.ru Igor Goryanin goryanin@gmail.com*

#### Specialty section:

*This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology*

Received: *27 November 2017* Accepted: *29 June 2018* Published: *25 July 2018*

#### Citation:

*Ivanov Y, Kozlov AF, Galiullin RA, Tatur VY, Ziborov VS, Ivanova ND, Pleshakova TO, Vesnin SG and Goryanin I (2018) Use of Microwave Radiometry to Monitor Thermal Denaturation of Albumin. Front. Physiol. 9:956. doi: 10.3389/fphys.2018.00956* This study monitored thermal denaturation of albumin using microwave radiometry. Brightness Temperature, derived from Microwave Emission (BTME) of an aqueous solution of bovine serum albumin (0.1 mM) was monitored in the microwave frequency range 3.8–4.2 GHz during denaturation of this protein at a temperature of 56◦C in a conical polypropylene cuvette. This method does not require fluorescent or radioactive labels. A microwave emission change of 1.5–2◦C in the BTME of aqueous albumin solution was found during its denaturation, without a corresponding change in the water temperature. Radio thermometry makes it possible to monitor protein denaturation kinetics, and the resulting rate constant for albumin denaturation was 0.2 ± 0.1 min−<sup>1</sup> , which corresponds well to rate constants obtained by other methods.

Keywords: water clusters, protein denaturation, rate constant, microwave emission, temperature

# INTRODUCTION

In medical studies, microwave radiometry (MWR) is based on the intensity of intrinsic electromagnetic radiation of the patient's tissues in the microwave range (Barrett and Myers, 1975; Carr, 1989; Leroy et al., 1998). MWR has found wide application in medical diagnostics (Carr, 1989; Leroy et al., 1998; Han et al., 2001; Stauffer et al., 2013, 2014; Livanos et al., 2018). To date, commercially available equipment has been developed to allow early diagnosis of cancers and other diseases (Vesnin, 2010) 1 . In the field of biochemical research, MWR is just starting to evolve (Ivanov et al., 2016a,b).

In this paper, the possibility of using MWR to study kinetics of protein denaturation is demonstrated. The study of thermal denaturation of proteins is important for understanding protein stability and for modeling their physicochemical properties. We used MWR to monitor the denaturation of albumin, the best characterized protein. Human serum albumin (HSA) is the most abundant protein in human blood plasma. Genes encoding HSA reside on chromosome 4 (Nishio and Dugaiczyk, 1996). Measurements of HSA are very important for diagnosis and treatment of many diseases (Anderson, 2000; Green et al., 2012). Since this protein is highly homologous to Bovine Serum Albumin (BSA) (Both are stabilized with 17 disulfide bonds and have the same

<sup>1</sup>Localized Cancer Tumor Detection Using Microwaves and Nanoparticles. Patent US15617440 registered (2017)-06-08, Priority date 2016-03-15.

isoelectric point.), the latter is used in modeling the properties of HSA (Gelamo et al., 2002; Vetri et al., 2007). BSA contains 583 amino acids (Steinhardt et al., 1971; Hirayama et al., 1990; Gelamo et al., 2002; Moriyama et al., 2008; Ahmad et al., 2011) and has a molecular weight of 66.4 kDa. Its secondary structure consists mainly of alpha-helices. Denaturation and protein aggregation are generally investigated using scanning calorimetry and optical methods (Vetri et al., 2007; Borzova et al., 2016). The BSA denaturation rate constant at a temperature of 60◦C is 0.5 min−<sup>1</sup> (Borzova et al., 2016). Thus, using BSA denaturation as an example, this paper demonstrates the use of MWR in biochemical studies. Existing biochemical methods, including those used to study albumins (Steinhardt et al., 1971) typically have a number of disadvantages. Optical methods, such as circular dichroism methods require use of optical spectral instrumentation, which is quite expensive and may be difficult to operate. The drawbacks of calorimetry are the complexity and duration of the analysis. It should be noted, the advantages of using new methods in addition to known methods are that they provide more information about the denaturation processes under investigation.

In this work, new applications of this technique for measuring kinetics of albumin denaturation are demonstrated by monitoring the change in the brightness temperature in the microwave range. In the long term, development of MWR and its adaptation for studying other biochemical processes, could potentially result in a new analytical platform in biomedical research. MWR is easy to operate, does not require labels, and as noted earlier, can be used to study radiative properties of biological systems.

# MATERIALS AND METHODS

De-ionized, ultrapure water (Milli-Q System, Millipore, USA) was used for all solutions. Lyophilized bovine serum albumin BSA [Sigma] was dissolved in deionized water at a concentration of 0.1 mM, close to the concentration of albumin in vertebrate blood (Fasano et al., 2005). pH of BSA solution was 7.00 ± 0.05.

#### MEASUREMENT OF RADIOMETRIC BRIGHTNESS IN THE MICROWAVE RANGE

We used a microwave radiometer RTM-01 RES (www.mmwr. co.uk) to measure electromagnetic radiation in the microwave range. The power of electromagnetic radiation in the frequency range 1f of the medium is proportional to the radiometric (brightness) temperature of the medium Trad (Weisblat, 2003). For the microwave range, the power of the microwave signal at the antenna output can be represented as:

$$P = kT\_{rad} \Delta f \, (1 - R) \,, \tag{1}$$

where 1f is the frequency band of the microwave radiometer, in GHz, and k is Boltzmann's constant, R is the reflection coefficient of antenna, P- is the microwave radiation power of the heated aqueous medium of the protein received by the antenna.

The antenna is designed in such a way as to minimize the reflection coefficient R. In addition, in the zero-balance radiometer used in the present studies, the measurement results are weakly dependent on the change in the reflection coefficient. Therefore, for balanced zero radiometers, R = 0 is usually assumed. In accordance with eq. (1), by measuring the power of the radiation of the medium in the microwave range, it is possible to obtain information on the brightness Trad(Brightness Temperature derived from Microwave Emission, BTME), which characterizes the emissivity of the water-protein medium. The measurement range of 3.4–4.2 GHz is selected based on the results of our early experiments (Ivanov et al., 2016a,b). So, in this range we observed the radiation of a solution of enzyme biomolecules in the course of their functioning. The functioning is accompanied by the restructuring of the hydrate shells of bio-molecules, so we assumed that the melting processes of the protein can also lead to a rearrangement of its hydration envelope, and, consequently, to changes in the spectral characteristics of the water-protein medium in this frequency range. In addition, the convenience of using this range is determined by the fact that it operates commercial radiometers. We measured the BTME of water and the BTME of an aqueous solution of BSA in the microwave range 3.4–4.2 GHz. The error in measuring the BTME was ±0.1◦C. To record radiation in the microwave range, a whip antenna connected to the microwave radiometer was used.

# PROCEDURE TO PREPARE SOLUTIONS FOR MEASUREMENT

Water and 0.1 mM aqueous protein solutions were heated using a thermostat (T-24, BIS, Russia) to a temperature of 56◦C, after which the sample was taken by pipette with a volume of 1 mL. Measurements of BTME were carried out after the solution was injected into a conical polypropylene measuring cuvette with a conical base diameter of 25 mm and a height of 17 mm. The cuvette was inserted into a thermostat Thermomixer Comfort (Eppendorf, Germany). Measurements were carried out for 25– 60 min. All measurements were repeated twice.

# RESULTS

The control series of the experiments included measurements of BTME (t) of pure water and BSA solution at a temperature of 23◦C, when BSA denaturation is not observed. In either case, there was no change in BTME (t) during the observation time of 22 min.

Radiothermometry was used to monitor kinetics of thermal denaturation of BSA in water at a temperature of T = 56◦C. We analyzed the BTME of an albumin solution in water, we have used the BTME of water as a control to exclude its possible effect on the results. During the observation time of 22 min, changes in the BTME albumin solution were observed (see **Figure 1**), while from 25 to 60 min of observation, water and the BSA solution were practically unchanged (therefore no data in this time range are presented). The BTME of water did not change

significantly with time (**Figure 1**). At the same time, the BTME of BSA changed significantly. During the initial period, the BTME of BSA was less than the BTME of water at ∼4 ◦C. During an observation time of about 25 min, BTME of BSA increased on the order of 1.5–2 ◦ C, (**Figure 1**). Error of experiment was of the order of 0.2◦C. The thermodynamic temperature of water and the albumin solution did not change during the experiment.

For better explanation of the experimental data, the effect of water was excluded, and a difference curve between the BTME of the protein solution and BTME of the water (control) was calculated. The difference between the BTME of BSA and that of water, 1BTME (t) = BTME (BSA)–BTME (water), is shown in **Figure 2**.

The 1BTME(t) of BSA decayed monotonically (**Figure 2**). This dependence was approximated by an empirical equation:

$$
\Delta BTME\left(t\right) = Ae^{-k\_{dm}t} + B,\tag{2}
$$

where A is the pre-exponential factor, kden is the denaturation rate constant, and t is the time of measurement.

The thermal denaturation curve of albumin (**Figure 2**) is well described by the exponential curve (2) with kden = 0.10 min−<sup>1</sup> , where A = 1.49, B = −3.77. Approximation error BTME (t) does not exceed 0.2◦C. Average value of kden = (0.2 ± 0.1) min−<sup>1</sup> . Note that our data showing significant changes in BTME of BSA at 56◦C is in agreement with experimental data derived from scanning calorimetry (Borzova et al., 2016), indicating that the BSA denaturation process should be observed in this temperature range. The denaturation rate constant at 60◦C, determined using scanning calorimetry, was <sup>k</sup>den ∼= 0.5 min−<sup>1</sup> (Borzova et al., 2016). Thus, the BTME (t) of BSA solution determined using MWR is close to that value. The reason for the BTME change in albumin solution during its thermal denaturation is as follows. First, it is known that while BSA denatures, the ratios

of ortho/para isomers of water change at the same time, as noted for chymotrypsin solution (Bunkin and Pershin, 2013). Second, as albumin denatures, it aggregates (Borzova et al., 2016), leading to a change in the ratio of ortho/para isomers of water (Pershin, 2013). The rotation spectrum of water isomers is in the GHz range (Bunkin and Pershin, 2013). These measurements are also made with a radiothermometer in the GHz band, so the denaturation of the protein leads to a change in the spectral properties of the solution within the sensitivity range of the instrument. Thus, MWR is well suited to monitoring protein denaturation.

MWR monitors the process of denaturation as a set of several processes described above, which allows refinement of models of protein denaturation. This method does not require any isotopes or fluorescent labels. It is simple and inexpensive compared with calorimetry and optical spectrometry and it allows quantum effects associated with the ratio of ortho/para isomers of water to be evaluated.

#### DISCUSSION

The work demonstrates that microwave radiometry allows direct monitoring of thermal denaturation of proteins, without the use of labels. This distinguishes it from complex methods using fluorescent labels or isotope tags. This method is much cheaper than circular dichroism or microcalorimetric methods. Within the microwave range, microwave radiometry makes it possible to monitor kinetics of protein denaturation associated with changes in these spectral properties of a medium containing a protein. These spectral properties are determined by the quantum-mechanical characteristics of this medium and provide information on the nature of the intermolecular interaction. Such characteristics include description the rotational transitions of water molecules to ortho- and para-states, the degree of disequilibrium of the aqueous medium in terms of its spin temperature, reflecting the complex heterostructure of water (Bunkin and Pershin, 2013). Moreover, some of the para-H2O molecules form hydrogen-bound complexes, for example, the structures of the hydrate shells of proteins (Pershin, 2013). The ratio of ortho/para isomers of pure water at constant temperature, as seen in **Figure 1**, does not change. However, in the presence of protein, the possible conversion of ortho/para isomers of water can be stimulated due to the change in the hydrate coat of the protein during its melting. The rotational transitions of hydroxy groups are in the GHz range, and the values of the detuning of the ortho/para states of the water molecules are also in this range (Pershin, 2013). In the same range, the spectral properties of the water-protein medium during the melting of the protein were measured. This allows us to suggest that changes in the brightness temperature of the medium measured in the range 3.4–4.2 GHz during the melting of the protein are associated with a change in the ratio ortho/para isomers an aqueous medium. Therefore, using this method, it is possible to observe quantum mechanical transitions of ortho/para isomers an aqueous medium in a protein solution during denaturation process. Accordingly, it is possible to obtain new information for the theoretical modeling of protein denaturation when the protein interacts with an aqueous microenvironment.

Microwave radiometry has a high potential to become a standard method for studying physico-chemical properties of protein media. This is because in addition to characterizing the absorbing properties of aqueous protein solutions, microwave radiometry makes it possible to determine the radiative properties of enzyme solutions (Ivanov et al., 2016a,b) and to

REFERENCES


Fasano, M., Curry, S., and Terreno, E. (2005). IUBMB. Life 57, 787–796.

Gelamo, E. L., Silva, C. H. T. P., Imasato, H., and Tabak, M. (2002). Interaction of bovine (BSA) and human (HSA) serum albumins with ionic surfactants: spectroscopy and modelling. Biochim. Biophys. Acta 1594, 84–99. doi: 10.1016/S0167-4838(01)00287-4

monitor cellular functions in the microwave range (Zynov'ev and Vesnin, 2009; Ivanov et al., 2016c).

# FUTURE APPLICATIONS IN BRAIN STROKE AND TRAUMA

The brightness temperature of water and protein solutions at physiological temperatures (Zynov'ev and Vesnin, 2009; Ivanov et al., 2016d, 2017; Vesnin et al., 2017), should permit applications this method in diagnostics of human and animal hemodynamics and for describing processes of water balance in living systems. In clinics, MWR has been successfully used for breast cancer (Vesnin et al., 2017), Among other clinical applications, MWR has long been used to monitor brain damage (Butrov et al., 2012) after stroke and other brain injuries. In the future, microwave radiometry could have clinical applications in stroke and trauma.

#### AUTHOR CONTRIBUTIONS

YI supervised the project. AK, RG, VT, VZ, and NI carried out the experiments. TP carried out the experiments and wrote the paper. SV provided equipment and consultancy. IG wrote the paper and provided consultancy.

### ACKNOWLEDGMENTS

This work was supported by the Program for Basic Research of State Academies of Sciences for 2013–2020.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Ivanov, Kozlov, Galiullin, Tatur, Ziborov, Ivanova, Pleshakova, Vesnin and Goryanin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Dynamic Cerebral Autoregulation Remains Stable During the Daytime (8 a.m. to 8 p.m.) in Healthy Adults

Wei-tong Guo<sup>1</sup> , Hongyin Ma<sup>1</sup> , Jia Liu<sup>2</sup> , Zhen-Ni Guo1,3 \* and Yi Yang1,3 \*

<sup>1</sup> Department of Neurology, The First Hospital of Jilin University, Changchun, China, <sup>2</sup> Institute of Advanced Computing and Digital Engineering, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, University Town of Shenzhen, Shenzhen, China, <sup>3</sup> Department of Neurology, Clinical Trial and Research Center for Stroke, The First Hospital of Jilin University, Changchun, China

Many functions of the human body possess a daily rhythm, disruptions of which often lead to disease. Dynamic cerebral autoregulation (dCA) stabilizes the cerebral blood flow to prompt normal neural function. However, whether dCA is stable across the day remains unknown. This study aimed to investigate the daily rhythm of dCA. Fifty-one healthy adults (38.294 ± 13.279 years, 40 females) were recruited and received six dCA measurements per individual that were conducted at predefined time points: 8:00, 9:00, 11:00, 14:00, 17:00, and 20:00. Although the blood pressure fluctuated significantly, there was no statistical difference in phase difference and gain (autoregulatory parameters) across the six time points. This study demonstrates that dCA remains stable during the interval from 8 a.m. to 8 p.m. and underscores the importance of stable dCA in maintaining cerebral blood flow and neural function.

#### Edited by:

Alexey Goltsov, Abertay University, United Kingdom

#### Reviewed by:

Giovanna Zoccoli, Università degli Studi di Bologna, Italy Evgeny Zherebtsov, Aston University, United Kingdom

#### \*Correspondence:

Zhen-Ni Guo zhen1ni2@163.com Yi Yang doctoryangyi@163.com; yang\_yi@jlu.edu.cn

#### Specialty section:

This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology

Received: 25 July 2018 Accepted: 30 October 2018 Published: 20 November 2018

#### Citation:

Guo W, Ma H, Liu J, Guo Z-N and Yang Y (2018) Dynamic Cerebral Autoregulation Remains Stable During the Daytime (8 a.m. to 8 p.m.) in Healthy Adults. Front. Physiol. 9:1642. doi: 10.3389/fphys.2018.01642 Keywords: dynamic cerebral autoregulation, blood pressure, cerebral blood flow, daily rhythm, cerebrovascular disease

# INTRODUCTION

Many functions of the human body have their own circadian rhythms and change over the course of a day. Dynamic cerebral autoregulation (dCA) is a process to maintain cerebral blood perfusion at an appropriate state via regulating cerebral vasculature (Budohoski et al., 2013). Cerebral autoregulation typically works with a mean arterial pressure (MAP) in the range of 60 and 150 mm Hg in normal conditions (Paulson et al., 1990). Previous studies have shown that patients with hemorrhagic stroke (Ma et al., 2016), ischemic stroke (Guo et al., 2014; Salinet et al., 2018), generalized anxiety (Guo et al., 2018), epilepsy (Lv et al., 2018), and Alzheimer's disease (Shekhar et al., 2017) have impaired dCA. Cerebral autoregulation can be affected many factors, such as CO<sup>2</sup> (Meng and Gelb, 2015), nitric oxide (Guo et al., 2016) and altitude level (Jansen et al., 2007). However, the daily rhythm of dCA and whether the time of measurement interferes with the evaluation of dCA are unclear. This study aimed to determine the rhythm of dCA during the interval from 8 a.m. to 8 p.m. in healthy adults.

### MATERIAL AND METHODS

#### Study Protocol

The study performed six dCA measurements at six predefined time points: 8:00, 9:00, 11:00, 14:00, 17:00, and 20:00. All the subjects had no history of chronic diseases or acute infections within the 2 weeks before beginning the study. Subjects with intracranial and/or extracranial major vascular Guo et al. Daily Rhythm of dCA

stenosis/occlusion diagnosed by a transcranial Doppler (EMS-9PB, Delica, China) and carotid ultrasound (IU22, Phillips, Andover, MA) were excluded. All the subjects were informed of the study protocol 1 day in advance. Subjects were instructed to refrain from caffeinated drinks and alcohol ingestion for 24 h before the examination (Claassen et al., 2016). All the subjects were non-smokers. The subjects were asked to avoid strenuous exercise, caffeine, and alcohol during the whole research process. The study was approved by the Ethics Committee of the First Hospital of Jilin University. Written consent was provided by all participants.

#### DCA Measurement

Measurements were performed in the specific examination room by a professional technician. The room was quiet and had a controlled temperature of 20◦C to 24◦C. The subjects were asked to take a relaxed supine position for 10 min. First, the technician measured the arterial blood pressure (ABP) at the brachial artery by an automatic blood pressure monitor (Omron 711, Japan). The continuous ABP was measured non-invasively using a servo-controlled plethysmograph (Finometer Model 1, FMS, Netherlands) at the middle finger. Simultaneously, two 2 mHz transcranial Doppler probes were placed over the temporal windows to monitor in real-time the bilateral middle cerebral arteries at a depth of 45–60 mm. the probes were fixed with a customized head frame to make sure cerebral blood flow velocity (CBFV) was continuously and stably measured. CBFV and continuous ABP were recorded simultaneously from each subject for 10 min. All data were recorded for further assessment and analysis.

## Data Analysis

Data of ABP and CBFV were acquired using MATLAB (MathWorks, Inc., US). The dynamic relationship between ABP and CBFV was analyzed by transfer function analysis (TFA) as follows (Claassen et al., 2016). A cross-correlation function between ABP and CBFV was used to align the data on in order to eliminate the possible time lags. An anti-alias filter, a third-order Butterworth low-pass filter, with cutoff frequency at 0.5 Hz was applied so as to down-sample the data to 1 Hz. Welch's method was employed to estimate the autospectrum of ABP, Sxx(f), and the cross-spectrum of ABP and CBFV, Sxy(f), in frequency domain by averaging the periodograms of the down-sampled ABP and CBFV with a 50% overlapped hamming window of 90 s. The transfer function, H(f), was then deviated as:

$$H\left(f\right) = \frac{\mathbb{S}\_{\text{xy}}(f)}{\mathbb{S}\_{\text{xx}}(f)}\tag{1}$$

Gain and phase difference (PD) can then be calculated from (1) by Equation (2) and (3), respectively, as:

$$\left|H\left(f\right)\right| = \sqrt{\left|\left|H\_R(f)\right|^2 + \left|H\_I(f)\right|^2\right|},\tag{2}$$

$$\theta = \tan^{-1} [H\_I(f) / H\_R(f)].\tag{3}$$

where R and I denote the real and imaginary parts of the transfer function, respectively. Phase difference (PD), gain, and coherence function within a 0.06–0.12 Hz frequency range were then derived from TFA to evaluate dCA. A low value of PD indicates that CBFV follows the changes of ABP passively, whereas a high value of PD suggests that CBFV is actively regulated against the fluctuations of ABP (van Beek et al., 2008). Due to the fact that TFA is a linear model-based method, signals with low coherence between ABP and CBFV (≤0.40) were excluded from further statistical analysis.

#### Statistical Analysis

The statistical analysis of the data was conducted using IBM SPSS Statistics 24.0 (Armonk, NY, United States), and a twotailed P-value < 0.05 was considered to be statistically significant. The distribution of data was assessed using a one-sample Kolmogorov—Smirnov test. Data are shown as mean ± standard deviation (SD) for normally distributed continuous variables. Repeated measurement analysis of variance was performed for comparing differences in observed values at different time points.

# RESULTS

The current study enrolled 51 healthy adults (38.294 ± 13.279 years, 40 females). DCA was measured in each subject for six times according to the preset procedure. Thus, the study included 306 records in total for dCA analysis. The coherence of all records was over 0.40. ABP and heart rate of serial measurements are presented in **Table 1**. ABP including systolic blood pressure (SBP), diastolic blood pressure (DBP), and mean arterial pressure (MAP) was statistically significant at different time points (P < 0.001, **Figure 1**, **Table 1**). Heart rate remained stable during the whole 12 h in a day (**Figure 1**, **Table 1**). The results of female or male independent analysis and whole subject analysis were consistent (**Table 1**).

### DCA

From 8:00 to 20:00, the PD was not significantly changed within 12 h (P = 0.233). Furthermore, the PD and gain showed no difference between left and right hemispheres at all the time points (P = 0.573, 0.388, 0.854, 0.263, 0.200, 0.665). Interestingly, the PD tended to be higher at 20:00 compared with the value at 8:00, but it was not statistically significant (**Table 1**, **Figure 1**). The gain did not significantly differ among all study time points (**Table 1**, **Figure 1**).

## DISCUSSION

The present study investigated the daily rhythm of dCA and the relationship between ABP and dCA. The major finding of the study was that dCA remains at a stable level during the daytime. Furthermore, dCA did not fluctuate following the changes of ABP.

Clinically, dCA measurement can correctly reflect most aspects of the autoregulatory response (Tiecks et al., 1995). Previous studies have shown that dCA was affected by the sympathetic nerve (Hamner et al., 2010), cholinergic nerve (Hamner et al., 2012), myogenic mechanisms (Tan et al., 2013), and metabolic control (Payne, 2018). Although these


SBP, systolic blood pressure; DBP, diastolic blood pressure; MAP, mean arterial pressure; HR, heart rate; cpm, counts per minutes.

TABLE 1 | Dynamic cerebral autoregulation

 parameter (phase difference, gain), mean arterial blood pressure, and heart rate and statistical results.

that we collected the data. (B) Statistical analysis of serial dCA, arterial blood pressure, and heart rate of serial measurements. The solid point represents means, and the whiskers denote standard deviation. \*P < 0.05 for comparison of different time points by repeated measurement analysis of variance. HR, indicates heart rate;

mechanisms have their own change rules, we showed that dCA is maintained at a stable level. The present study may indicate that the daily rhythm of dCA is independent of the change rule of

cpm, counts per minutes; SBP, systolic blood pressure; DBP, diastolic blood pressure; MAP, mean arterial pressure.

one single mechanism. Cerebral autoregulation functions to stabilize cerebral blood flow and metabolism while the ABP changes (Lassen, 1959). When dCA was impaired, the cerebral blood flow would change with fluctuations in blood pressure, leading to cerebrovascular disease. ABP fluctuation has been found to be a potent risk factor for cerebrovascular disease (Yano and Kario, 2012). Due to the morning blood pressure surge, most strokes occur in the morning hours (Marler et al., 1989; Yano and Kario, 2012). Similar to previous studies (Atkinson et al., 2010; Douma and Gumz, 2018), the present study reproduced the finding of ABP fluctuations. However, we found that dCA was independent of ABP fluctuations during the daytime.

Our analysis indicated that healthy adults have stable and functional dCA to maintain brain perfusion, similar to the findings of a previous study (Brys et al., 2003). Chi et al have demonstrated that dCA was interchangeable and effective by assessing from the internal carotid artery compared with that from the middle cerebral artery (Chi et al., 2017). Also, dCA assessed by using a 5 min recording was identical with that using a 10 min recording in the clinical application (Chi et al., 2018). Despite the stability of different recording time and different measuring locations, we have discovered the stable rhythm of dCA during the 8:00–20:00 interval. And, conversely to the abundant literature showing differences in cerebral blood flow between sexes (Barnes, 2017) and the influence of sex on circadian rhythms (Yan and Silver, 2016), the present study has shown the same results for both sexes.

In the literature, dCA was reported to be reduced at 6 a.m.−8 a.m. (Ainslie et al., 2007), which is different from the present study. We speculate that the differences might arise from three reasons: (1) different detection time of cerebral autoregulation: Ainslie et al performed the morning dCA measurement in the early morning (6 a.m.−8 a.m.), while our study performed the morning dCA measurement at 8 AM. Because we don't know the differences between "early morning" and "morning" of the dCA, we cannot exclude the difference caused by the detection time points; (2) the sex of the subjects: Ainslie et al recruited all male subjects, while the present study recruited both male and female subjects; (3) the age of the subjects; Ainslie et al recruited young subjects with a mean age of 25. However, our study recruited middle-aged subjects with a mean age of 38.

The present study shows that dCA is a reliable indicator and remains at a stable level during the daytime and was not affected by fluctuations of ABP. Thus, given its stability and reliability, it possible to assist in diagnosing and assessing the cerebrovascular function of cerebrovascular diseases, generalized anxiety, epilepsy and Alzheimer's disease. Simultaneously, the dCA also can be used to evaluate the treatment effect of the above disease.

However, there are several limitations to the study. The current study lacks the nocturnal and early morning rhythm of dCA. The main reason is that monitoring dCA requires the subjects to stay awake. If night monitoring is carried out, subjects may not sleep well, which may affect daytime results. Furthermore, due to the nature of the observational study, we did not collect the blood samples at different time points to test the changes of biomarkers. These questions should be pursued in future studies.

#### REFERENCES


In summary, the rhythm of dCA keep was steady during the 8 a.m.−8 a.m. interval in healthy adults, and it is not influenced by the fluctuations of ABP. For evaluating dCA, one random measurement of dCA is reliable.

#### AUTHOR CONTRIBUTIONS

WG design the study, acquired the data, analyzed and interpreted of data, drafted the manuscript. Z-NG analyzed and interpreted of data and revised the manuscript. JL analyzed and interpreted of data and obtained funding. HM acquired the data and revised the manuscript. YY studied concept and designed the study, critical revised the manuscript, supervised the study and obtained funding.

#### FUNDING

This article was supported by the National Key R&D Program of China (2016YFC1301600), JLUSTIRT (2017TD-12), and Jilin provincial development and reform commission (2017C018) to YY, and National Natural Science Foundation of China (Grant NO. 81871447) to JL.

#### ACKNOWLEDGMENTS

We thanked all the volunteers for their contributions to the study.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Guo, Ma, Liu, Guo and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Multimodal Optical Diagnostics of the Microhaemodynamics in Upper and Lower Limbs

Angelina I. Zherebtsova<sup>1</sup> , Viktor V. Dremin<sup>1</sup> , Irina N. Makovik<sup>1</sup> , Evgeny A. Zherebtsov1,2,3 \*, Andrey V. Dunaev<sup>1</sup> , Alexey Goltsov<sup>4</sup> \*, Sergei G. Sokolovski3,5 and Edik U. Rafailov3,5

<sup>1</sup> Research and Development Center of Biomedical Photonics, Orel State University, Oryol, Russia, <sup>2</sup> Optoelectronics and Measurement Techniques Unit, University of Oulu, Oulu, Finland, <sup>3</sup> Optoelectronics and Biomedical Photonics Group, Aston Institute of Photonic Technologies, School of Engineering and Applied Science, Aston University, Birmingham, United Kingdom, <sup>4</sup> School of Applied Sciences, Abertay University, Dundee, United Kingdom, <sup>5</sup> International Center of Critical Technologies in Medicine, Saratov State University, Saratov, Russia

#### Edited by:

Antonio Colantuoni, University of Naples Federico II, Italy

#### Reviewed by:

Claudia Penna, University of Turin, Italy Dominga Lapi, University of Pisa, Italy

#### \*Correspondence:

Evgeny A. Zherebtsov evgenii.zherebtsov@oulu.fi Alexey Goltsov a.goltsov@abertay.ac.uk; alexey.goltsov@gmail.com

#### Specialty section:

This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology

Received: 01 November 2018 Accepted: 27 March 2019 Published: 16 April 2019

#### Citation:

Zherebtsova AI, Dremin VV, Makovik IN, Zherebtsov EA, Dunaev AV, Goltsov A, Sokolovski SG and Rafailov EU (2019) Multimodal Optical Diagnostics of the Microhaemodynamics in Upper and Lower Limbs. Front. Physiol. 10:416. doi: 10.3389/fphys.2019.00416 The introduction of optical non-invasive diagnostic methods into clinical practice can substantially advance in the detection of early microcirculatory disorders in patients with different diseases. This paper is devoted to the development and application of the optical non-invasive diagnostic approach for the detection and evaluation of the severity of microcirculatory and metabolic disorders in rheumatic diseases and diabetes mellitus. The proposed methods include the joint use of laser Doppler flowmetry, absorption spectroscopy and fluorescence spectroscopy in combination with functional tests. This technique showed the high diagnostic importance for the detection of disturbances in peripheral microhaemodynamics. These methods have been successfully tested as additional diagnostic techniques in the field of rheumatology and endocrinology. The sensitivity and specificity of the proposed diagnostic procedures have been evaluated.

Keywords: laser Doppler flowmetry, tissue reflectance oximetry, pulse oximetry, fluorescence spectroscopy, skin thermometry, rheumatic diseases, diabetes mellitus

# INTRODUCTION

Blood microcirculation plays an important role in the transport of nutrients, oxygen, hormones and the release of metabolic products. In the last decade, there has been a steady increase in the interest of researchers in the problems of microcirculatory disorders in patients with rheumatological and endocrinological diseases. This is due to the significant role of microcirculation in the pathogenesis of such diseases (Avouac et al., 2011; Sena et al., 2013; Fuchs et al., 2017).

The medical, social and economic significance of rheumatic and endocrinological diseases is determined primarily by their prevalence, as well as the development of complications in the majority of patients, which significantly reduce the duration and quality of life, up to disability. According to European League Against Rheumatism (EULAR) statistics, more than 22% of the population already has or had previously rheumatic diseases (Smolen et al., 2017). At the same time, according to the International Diabetes Federation (IDF), by 2017 the prevalence of diabetes mellitus in the world has reached 8.8% among the adult population (Karuranga et al., 2017).

Rheumatic diseases lead to morphological disturbances of the microcirculatory bed, which include the rarefaction of the capillary network, the asymmetry of the capillaries, and the

**88**

appearance of megacapillaries (up to 50 µm in diameter) (Gutierrez et al., 2010). Chronic hyperglycemia and insulin resistance in diabetes mellitus cause increased vascular permeability, disruption of autoregulation of blood flow and vascular tone (Avogaro et al., 2011; Kolluru et al., 2012), leading to the structural and functional changes in capillaries and arterioles (Schramm et al., 2006).

In the absence of adequate therapy, peripheral blood flow disorders in limbs can lead to painful ulceration, gangrene, and the development of cardiovascular diseases. Previously published studies have shown that the risk of developing cardiovascular disease in patients with rheumatic diseases is comparable to that in patients with type 2 diabetes (Peters et al., 2010; Kerekes et al., 2012).

The earliest, usually reversible manifestation of these diseases is the development of microcirculatory dysfunction (Vanhoutte et al., 2015; Fuchs et al., 2017) due to endothelial damage, excessive expression of certain adhesion molecules and other factors. At present, several invasive and noninvasive methods are used in the clinical practice of endocrinologists and rheumatologists, including colorimetric duplex dopplerography, angiography, computed tomography, magnetic resonance angiography, transcutaneous oximetry (TcpO2), etc. (Daly Susan and Leahy Martin, 2012).

The methods of optical diagnostics are promising for the study of early microcirculatory disorders in patients with rheumatic diseases and diabetes. These methods have several advantages: painless procedures, quick results, lack of expensive consumables, minimal impact on the object and its properties. This article presents an overview on recent advances in optical non-invasive diagnostics of the peripheral hemodynamics of the upper and lower limbs in rheumatological and endocrinological profile patients. Among these are laser Doppler flowmetry (LDF), tissue reflectance oximetry (TRO), pulse oximetry (PO) and fluorescent spectroscopy (FS).

The laser Doppler flowmetry method allows for investigating the blood flow in the microcirculatory bed in vivo. The method is based on probing the tissue with laser radiation and analyzing back reflected from the tissue radiation partially scattered from moving red blood cells (Bonner and Nossal, 1990; Krupatkin and Sidorov, 2013).

The tissue reflectance oximetry method provides information about the tissue oxygen saturation (StO2) of the examined biological tissue microcirculation and allows for calculation of the relative blood volume (Vb) in the surface layers of soft tissues (skin, mucous membranes). This technology is based on the ability of oxygenated and deoxygenated hemoglobin to absorb light in the red and near infrared range (Casavola et al., 1999; Wallace et al., 2009).

Mechanisms that cause abnormal microcirculation in pathological processes can be quantified by use of wavelet transform for the records of LDF and TRO signals. The amplitude-frequency analysis of oscillations in skin blood flow allows distinguishing five frequency ranges corresponding to metabolic (endothelial) (Kvandal et al., 2003), neurogenic, myogenic (Söderström et al., 2003; Krupatkin, 2006), respiratory and cardiac (Stefanovska et al., 1999; Krupatkin, 2008) activity.

Optical diagnostics enables not only evaluation of blood flow parameters but registration of associated biochemical changes in the living tissue under study. The fluorescence spectroscopy is one of such methods allowing for registration metabolic changes in vivo. As an example, the method has promises to be applied in diagnostics of the diabetes mellitus complications. An important indicator of the viability of tissues is the mitochondrial function. By the parameters of the respiratory chain, one can speak of normal or pathological activity of the cells, diagnose the state of tissue ischemia. One of the estimates of the mitochondrial function is the ratio of coenzymes NADH and FAD, which can be calculated from the intensity of their fluorescence (Bartolomé and Abramov, 2015). In addition, in recent years, it has been found that the long-term effects of pathogenic factors such as hyperglycemia and oxidative stress in diabetes mellitus (Tan et al., 2002) and chronic inflammatory diseases (Murdaca et al., 2012) can lead to increased glycation of proteins and accumulation of advanced glycation end products (AGEs), which affect the properties of collagen and other structural proteins of the capillary membrane and skin (Gkogkolou and Böhm, 2012). These changes can be quantified by the intrinsic fluorescence of the pentosidine residues formed during glycation of collagen (Sell et al., 1991; Takahashi et al., 1995).

To increase the reliability of the results obtained by the optical spectroscopy methods, an assessment of the measured parameters changes is carried out with functional tests. Occlusion, heat, cold, respiratory, orthostatic, electrostimulation et al. tests are currently used as an act of provocation. Functional tests allow for revealing latent hemodynamic disorders and adaptive reserves of the microcirculation system.

The recent studies show that multiparametric approach in the form of joint measurements by various optical methods supplemented by functional tests can give complex information about the functional state of the microcirculation (Goltsov et al., 2017; Zherebtsov et al., 2017; Zherebtsova et al., 2017).

The aim of this work is to generalize the previously obtained results (Dremin et al., 2017; Mizeva et al., 2017, 2018; Zherebtsov et al., 2017; Zherebtsova et al., 2017) and to assess the possibilities of using optical non-invasive methods in studying the microcirculatory bed of patients with disorders in peripheral microhemodynamics.

The article provides a review of three methods of combined use of optical non-invasive technologies for the diagnostics of blood microcirculation and metabolic disorders. The main distinctive feature of all described methods is the use of functional tests, which greatly enhanced the diagnostic capabilities of the methods (Herrick and Clark, 1998).

The first two experimental studies demonstrate the possibility of using optical non-invasive methods and cutaneous temperature for detection violations of peripheral blood flow of the upper limbs in rheumatological profile patients. The cold water exposure, occlusion tests (OT) and cold pressor test (CPT) were selected as a provocative actions. Occlusion test is most commonly used functional test to investigate and assess microvascular function according to post occlusive reactive hyperemia (PORH) response (Roustit and Cracowski, 2012). A cold test is also often used for diagnostic purposes due to

vasospasm is characterized by sporadic manifestation and influenced by trigger factors like cold exposure, emotional stress, physical exercise etc. (Terada et al., 2007; Ye and Griffin, 2016). The deliberate provocation of vasospasm allows one to increase the sensitivity of the diagnosis and to assess the severity of the pathological process.

The third experimental study estimates the potential of synchronous registering the blood flow parameters and the fluorescence of intrinsic tissue fluorophore with the purpose of diagnosing the stages of complications of lower limbs in diabetes mellitus patients. The local heating stimulation was chosen as test action on the blood microcirculation system. The local heating test allows the assessment of the local regulatory mechanisms of blood flow.

### SIMULTANEOUS MEASUREMENTS OF THE BLOOD PERFUSION AND SKIN TEMPERATURE FOR FUNCTIONAL DIAGNOSTICS OF INTRADERMAL FINGER VESSELS

#### Materials and Methods

The aim of this study was evaluation of the combined use of the laser Doppler flowmetry and skin thermometry methods during the occlusion test to distinguish vasospastic disorders in hands (Zherebtsov et al., 2017). The experimental studies involved 27 healthy volunteers (HV) (mean age 23 ± 5 years) and 41 patients with rheumatic diseases (PRD) (mean age 56 ± 12 years) from the Rheumatology Department of the Orel Regional Clinical Hospital (Oryol, Russia). Subjects in the HV group did not have diagnosed diseases that are accompanied by secondary vasospastic syndrome, and did not have signs and symptoms which are indicative for primary vasospasm or predisposition to vasospasm [e.g., cold hands and feet, low blood pressure (Flammer et al., 2001)]. Presence of vasospasm of each tested subject in the PRD group was confirmed by the attending physician.

It is well known that microcirculatory abnormalities contribute to the pathogenesis and pathophysiology of numerous rheumatic diseases (Murray et al., 2004). Microvascular changes taking place in the early stages of such diseases may not manifest itself. In many cases, vascular problems can remain hidden. Such a high risk factor as aging can significantly contribute to the vascular state of the volunteers (Abularrage et al., 2005). For this reason, before carrying out measurements in the group of patients with diagnosed rheumatic diseases, we have tried to find a control group of people in their 50s without clear signs of vascular diseases. Nevertheless, due to the lack of confidence on absence of vascular disorders among the aged volunteers, the younger people who are more likely have no microcirculatory disorders were intentionally included in the control group.

The laser analyser of blood microcirculation "LAKK-02" (SPE "LAZMA" Ltd., Russia) and a custom developed multi-channel thermometry device for low inertia thermometry were used for experimental measurements. The measurements of cutaneous temperature and the index of microcirculation were performed on the distal phalanx of the third finger of the right hand. The specially designed attachment was used for longitudinal arrangement of the LDF fiber probe. The experimental setup and design features and location of the proposed attachment is presented in **Figure 1**.

The stage of heating of the hand at 42◦C was applied for secure equal initial conditions of the experiment. This temperature ensures complete dilatation of vessels and the rate of blood flow returns to a normal level, even in patients with vasospastic disorders. The temperature of the main part of experiment was 25◦C. It was selected from considerations of comfort for

patients and maximal potentially possible dynamic range of the skin postocclusion temperature response. Occlusion test was performed on the upper arm using a sphygmomanometer air cuff with a pressure of 200–220 mmHg for 3 min. The experiment was conducted according to the study protocol is described in **Figure 2**.

All measurements were carried out mainly in the morning, 2–3 h after meal, in a state of mental and physical rest. The patient sat in such a way that the forearm of his right hand was 20 cm below the level of the heart (**Figure 1A**). The total duration of the experiment was not more than 40 min. Representative record of the cutaneous blood perfusion and temperature of the conditionally healthy volunteer is shown in **Figures 3A,B**, respectively.

#### Data Analysis and Diagnostic Criteria

Based on the results of the study, it was revealed several characteristic types of the blood microcirculation system response to temperature and occlusion effects (Zherebtsova et al., 2015).

The first type of response is characterized by the absence of the spasm, that reflects in temperature and perfusion increased in post-occlusion period.

In the case of the second type of the response, the initial parameters of the thermoregulatory system are such that at a moderate cooling of fingers the system can rapidly switch to state with low temperature. However, the effect of vasodilation induced by the occlusion is higher than the vasoconstriction induced by cooling. Thus, the mechanism of endothelial function for patients with the second type of response remains active.

In the case of the third type of response, vasospasm tendency remains much high (a tendency to stay in the stationary point of low temperature) and the vasodilatory effect of the occlusion is minimal. Therefore, the vasodilation does not appear in this case against a background of the vasoconstriction of vessels for the present temperature level.

The analysis of the obtained results showed that for the reliable identification of the functional state of blood microcirculation in fingers it is advisable to use two parameters based on the type of response of the cutaneous temperature and the skin blood flow to the arterial occlusion at a low temperature.

For the assessment of the vasodilation effect on cutaneous temperature, it was proposed to estimate the ratio of the difference between the maximum temperature after occlusion and the temperature of cold water to the difference between the minimum cutaneous temperature during occlusion and the temperature of cold water. To account for the effect of biological tissue heat capacity, the obtained value was normalized to the volume of the distal phalanx of the finger. Thus, the index of temperature response (ITR, arb. un.), based on skin temperature measurements, was calculated by the following formula:

$$ITR = \frac{1}{V} \frac{(T\_{\rm PO} - T\_{\rm CW})}{(T\_{\rm O} - T\_{\rm CW})} \tag{1}$$

where T<sup>O</sup> – the minimum temperature of the biological tissue during occlusion period, ◦C; TPO – the maximum temperature of the biological tissue during occlusion period, ◦C; TCW – the temperature of cold water in the heat-insulated water bath; V – the volume of the distal phalanx of the test finger, calculated as the volume of semi-ellipsoid, cm<sup>3</sup> .

ITR allows assessing the extent of post-occlusion microvessels vasodilatation at artificially created low ambient temperature. Vasospasm can be described as a state when thermoregulatory system tends to the stable stationary point that is close to ambient temperature ("stationary point with low temperature"). Whereas in case of normality the thermoregulation is in a stationary point which is characterized by a temperature above the ambient temperature ("stationary point with high temperature"). Formation of the feedforward trigger system of switching between the stationary states is caused by the positive non-linear relationship between temperature and blood perfusion (Nakamura, 2011). Decrease in the temperature of the biological tissue leads to decrease in the blood perfusion, which in turn contributes to the further temperature reduction. The effect of vasoconstriction due to temperature falling can be compensated by vasodilation effect of occlusion (Dezfulian et al., 2017), which appears in normal conditions and characterizes the functional state of vascular endothelium in norm (Meredith et al., 1996). Decrease of index of the temperature response was observed in the PRD group.

However, the use of parameters from only cutaneous thermometry does not allow the relevant separation of the norm and presence of vasospastic disorders. In this connection, we suggest to use the composite diagnostic criteria, which includes both LDF and cutaneous thermometry parameters.

For the assessment of adaptation reserves in the blood microcirculation system and for the evaluation the process of skin blood flow restoration after occlusion, the most commonly used value is the percentage ratio between the maximum skin blood perfusion after occlusion and the average skin blood perfusion level in basal conditions. In these experiments, we used

the parameter called blood flow reserve (BFR, %), calculated according to the formula:

$$BFR = \frac{Im\_{\text{max}}}{Im\_{\text{base}}} 100\% \tag{2}$$

where Immax – average index of blood microcirculation in the first 60 s after occlusion, PU; Imbase – average index of blood microcirculation during 60 s before occlusion, PU. Averaging of skin blood flow during 60 s is used to exclude the influence of motion artifacts on the results of calculations.

Post occlusive reactive hyperemia is used to investigate and assess microvascular function (endothelial function). Therefore, the BFR parameter characterizing the reaction of microvessels to arterial occlusion can serve as an effective tool for assessing microcirculatory disorders in rheumatic diseases.

In order to synthesize the decision rule for diagnosis, it is necessary to find a discriminant function that separates classes of the norm and the presence of vasospastic disorders using the ITR and BFR parameters.

Experimental studies have shown that PRD have lower values of the blood flow reserve, as well as a reduced index of temperature response after occlusion. There is a statistically significant difference of BFR and ITR parameters between values calculated for the HV and PRD. Comparison of these parameters between studied groups is shown in **Figure 4**.

The values of BFR and ITR have a two-dimensional normal distribution. The training sample size n = n<sup>1</sup> + n<sup>2</sup> = 68 is more than 20 times larger than the number of variables in the vector of informative parameters m = 2.

All these conditions enable us to apply the mathematical method of linear discriminant analysis (Demir and Ozmehmet, 2005). Using this method a discriminant function was defined in a linear form

$$F(BFR, ITR) = 0.022 \cdot BFR + 1.61 \cdot ITR \tag{3}$$

that allowed synthesis of the desired decision rule:

$$\begin{cases} \text{healthy} & \text{if } F(\text{BFR}, ITR) > 3.7, \\ \text{vasospin disorder} & \text{if } F(\text{BFR}, ITR) \le 3.7. \end{cases} \tag{4}$$

Substituting the experimental values of the parameters BFR and ITR into eq. (3), one can determine the presence or absence of vasospastic disorders of studied subject using decision rule (4).

**Figure 5A** shows the scatter plot of parameters ITR and BFR with the discriminant function (3). Points of the given graph correspond to a combination of experimental values BFR and ITR for examined persons, and discriminant linear function (3)

divides the feature space into two half-planes. According to the decision rule (4), the area above the discriminant straight line corresponds to the absence of vasospastic disorders in the fingers, the area below – to the presence of vasospastic disorders.

Leave-one-out cross-validation method was used for verification of the decision rule (Rao et al., 2007). Verification procedure reveals the probability of false-negative diagnosis result at the level of 0.13 (sensitivity 0.87), false-positive at 0.26 (specificity 0.74). The receiver operating characteristic curve (ROC-curve) for the suggested method of diagnosis is shown in **Figure 5B**. It indicates the quality of the suggested method of binary classification. The curve shows the sensitivity to specificity ratio at various values of the threshold between the domains of normal state of health and angiospastic disorders. The area under the curve is 0.88, which is indicative of a high efficiency of the classifier.

The proposed decision rule (4) allows us to offer a diagnostics method for the functional state of peripheral vessels, which can be characterized as reserve possibilities of blood flow using the LDF method, and reactivity of peripheral vessels located at a greater depth using the cutaneous thermometry method.

# Experimental Results and Discussion

The new diagnostic procedure using methods of laser Doppler flowmetry and thermometry during combined provocative factors of cold and ischemia was proposed in the study. It was suggested to use the measurements of blood microcirculation and skin temperature during the occlusion test in the thermally stabilized environment as a diagnostic approach for identification dysfunctions of the peripheral blood vasculature in PRD.

Using the experimental values of the parameters of LDF- and thermograms, a simple model for the classification of the presence or absence of vasospastic disorders in the fingers was proposed. The resulting classification model showed good results of sensitivity (0.87) and specificity (0.74).

The proposed approach of using a combination of several diagnostic technologies makes it possible to conduct a comprehensive assessment of the functional state of microcirculation vessels and larger vessels of the fingers. The combined effect of low temperature and temporary circulatory arrest provokes a manifestation of the vasospastic state. Since the vasospasm appears sporadically, this approach allows reducing the level of false–negative diagnostic results.

The results from this conducted research can be used in the development of multi-functional non-invasive diagnostic systems for the diagnosis and prevention of diseases associated with changes in the functional state of peripheral vessels.

## COLD PRESSOR TEST IN DETECTION OF DISORDERS IN THE MICROCIRCULATORY BED OF UPPER LIMBS

#### Materials and Methods

To evaluate the combined use functionality of the laser Doppler flowmetry (LDF), tissue reflectance oximetry (TRO), pulse oximetry (PO) methods and cold pressor test (CPT), experimental studies involved 32 HV (mean age 22 ± 2 years) and 60 PRD (mean age 55 ± 14 years) from the Rheumatology Department of the Orel Regional Clinical Hospital (Oryol, Russia).

The group of patients includes individuals primarily with the rheumatoid arthritis and systemic lupus erythematosus.

As described previously, disorders of upper limb microcirculatory bed are most commonly found as one of the forms of rheumatological profile disease pathologies. These diseases are more common in the elderly. It is thus necessary to clearly differentiate between a healthy state and one with microcirculatory bed disorders. A group of healthy young volunteers was recruited as a control to ensure an "extreme" state of good health, as they would present the lowest chance of exhibiting any undesired physiological conditions.

The experiment was conducted according to the study protocol described in **Figure 6**.

All measurements were performed in conditions of physical and mental rest 2 h after a meal. Volunteers also underwent a preliminary adaptation to room temperature 24–25◦C for 15–20 min in a sitting position, with the right arm on the table at heart level. The adaptation of volunteers to standard room temperature and abidance of study protocol during all measurements reduces the influence of different factors on results of diagnosis.

Experimental systems "LAKK-OP" and "LAKK-M" (SPE "LAZMA" Ltd., Russia) were applied for the measurement LDF-, TRO-, and PO-signals. These diagnostic devices utilized identical measurement channels. Temperature of water in a container during the cold exposure was controlled by a contactless digital thermometer (Sensitec NB401, Netherlands). Schemes of experimental installation during the measurement, during the cold exposure and schematic presentation of the LDF- and TRO-probe and pulse oximetry sensor positioning on a finger are showed in **Figure 7**.

The used TRO channel calculates the tissue oxygen saturation, which is defined as the percentage composition of oxyhaemoglobin in the sum of major fractions of hemoglobin in a tissue volume (Dunaev et al., 2014). The instrument implements the computational model based on the modified diffuse approximation of the light transfer equation (Spott et al., 1997) and utilizes the differences in the spectral characteristics of oxygenated and deoxygenated hemoglobin. The TRO measurements allow for non-invasive monitoring of microhaemodynamics and transport and utilization of oxygen within the blood microcirculation system.

Basic microcirculatory bed parameters were registered during experimental studies, providing a vast array of information. Index of blood microcirculation (Im), tissue oxygen saturation (StO2) and arterial oxygen saturation (SàO2) were registered by proposed methods.

Wavelet analysis of the registered LDF- and TRO-signals was carried out to evaluate the regulatory mechanisms. The complex Morlet wavelet was used as the analyzing wavelet (Frick et al., 2015). The calculation of wavelet coefficients for the frequency range from 0.01 to 2 Hz was performed with a logarithmic partitioning into 50 frequency sub-bands. Global wavelet power spectra were also calculated for both study groups in each BT. Maximum amplitude of peripheral blood flow oscillations in one of the frequency bands [endothelial (Ae), neurogenic (An), myogenic (Am), respiratory (Ar) and cardiac (Ac)] calculated from the wavelet analysis of LDF- and TRO-signals (Goltsov et al., 2017). Based on measured parameters and results of wavelet analysis of the registered LDF- and TRO-signals using proposed approach (Makovik et al., 2017) complex parameters were calculated.

Statistical analysis of the measured and calculated parameters was performed using non-parametric criteria: the Mann–Whitney test for comparing values between groups and the Wilcoxon test for comparing values within a single group.

#### Data Analysis and Diagnostic Criteria

Results of experimental studies showed that CPT evokes different reaction from the microcirculatory bed in each group. In particular, partial or complete recovery of the blood flow parameters after the CPT in some subjects was observed. As it is shown in **Figure 8** higher perfusion in basal state in PRD

group was observed, herewith cooling provokes the more evident response to changes of temperature in HV (Mizeva et al., 2017).

Interestingly, that analysis of averaged wavelet spectra in BT1 showed that larger energy in high frequency pulsations (above the frequency of 0.24 Hz) of the blood flow were observed in PRD than in ÍV, herewith this energy saved and after the cooling. The difference between groups in the low-frequency part of the spectra does not observe. After 20 min total restore of blood flow as far as spectral composition was observed.

As was previously stated in the works (Gutierrez et al., 2010; Yildiz, 2010) this reaction is associated with a weak damping capacity of the vascular bed due to a decrease in the elasticity of the vascular wall and increase its stiffness, and also because of morphological disturbances arising during RD (formation of megacapillaries, thinning of the capillary network).

Based on the obtained results and the differences between HV and PRD, the values of the perfusion and the maximum amplitude of LDF oscillations during CPT was used to the

synthesis of the decision rule for diagnose microcirculatory disorders in RD. These parameters satisfy the principles of statistical independence, as well as the significance of the differences of their values, calculated for the PRD and HV. A discriminant function of the values of perfusion Im2 and the maximum amplitude of blood flow oscillation in cardiac frequency band Ac2 for BT2 (measurement immediately after the cold exposure) was defined in the following linear form

$$F(I\_{\rm m2}, A\_{\rm c2}) = 0.12 \cdot I\_{\rm m2} + 1.93 \cdot A\_{\rm c2} - 3.25 \tag{5}$$

**Figure 9A** shows the scatter plot of parameters Im2 and MT<sup>2</sup> with the discriminant function (5). Points of the given graph correspond to a combination of experimental values Im2 and MT<sup>2</sup> for HV and PRD, and discriminant linear function (5) divides the feature space into two half-planes. The area above the discriminant straight line corresponds to the absence of microcirculatory disorders in the fingers, the area below – to the presence of microcirculatory disorders.

**Figure 9B** shows the ROC curve calculated for the obtained discriminant function. Area under curve (AUC) was used to compare the quality of different classifying rules. For the synthesized decision rule AUC equals 0.92. The results demonstrate that the perfusion and amplitude of the pulse wave can act as independent markers for microcirculatory disorders in RD.

Additional calculated parameters have been evaluated as it described in Makovik et al. (2017). Analysis of the parameters in each group during experimental study identified differences in myogenic tone (MT) and rate of oxygen consumption (OC) in BT3 (measurement 15 min after the cold exposure) (**Figure 10**), namely higher level of MT and lower level of OC in PRD than in HV.

Differences in OC level in PRD and HV (**Figure 10B**) indicate the possible violations of the microvascular bed surface of the smallest arterioles and capillaries (Makovik et al., 2017). In PRD a decrease of OC with an increase of MT compared with the values for HV are observed (**Figure 10A**). Such result can be interpreted as a sign of reduction of oxygen diffused through the vessel walls. These processes can lead to hypoxia, edema of tissues and the appearance of necrobiotic processes at untimely diagnostics and absence of treatment.

An analysis of the possible causes of these pathological changes revealed their association with an increase in myogenic tone, as well as its combination with venous stasis. So MT and OC can be used as an additional diagnostic criteria for detection of complications associated with microvascular disturbances and their possible causes.

#### Experimental Results and Discussion

Experimental studies have shown that the combined use of optical non-invasive technologies (laser Doppler flowmetry and tissue reflectance oximetry) in combination with cold pressor test has a huge diagnostic potential. The described method allows diagnosing the presence of vasospastic disorders in rheumatic diseases and identifying the possible cause of the pathological condition.

The proposed approach is based on the use of LDF and TRO methods before and after the exposure and subsequent wavelet transform of the signals. The level of the skin blood flow and its spectral properties undergo changes due to violations in the microcirculation system regulation. The mean level of the local microcirculatory blood perfusion and the amplitude of the pulse oscillations of blood flow immediately after CPT are two proposed diagnostic criteria, which were chosen for the synthesis of the decision rule. The resulting classification model provides excellent results of sensitivity (0.92) and specificity (0.97) in diagnosing microvascular disorders in rheumatic diseases (see **Figure 9B**).

In case of detection of microcirculatory disorders, the second stage of the proposed diagnostic method was implemented. It consists in identification of the associated complications and their possible causes with the use of additional diagnostic criteria: complex parameters of hemodynamics (myogenic tone, MT) and tissue respiration (rate of oxygen consumption, OC). The parameters are calculated based on the measured LDF and

the classification method.

TRO parameters before and after CPT. The combination of parameter values allows assessing the presence of complications and identifying myogenic or myogenic-congestion cause of their occurrence.

The clinical trail of the method showed a considerable promise in clinical application for the functional evaluation of microcirculatory blood flow regulation and has good possibilities in the diagnosis of microvascular disorders associated with rheumatic diseases.

# MULTIMODAL OPTICAL MEASUREMENTS UNDER LOCAL HEATING TEST TO REVEAL CHANGES IN LOWER LIMBS IN PATIENTS WITH DIABETES MELLITUS

#### Materials and Methods

The aim of this study was to estimate experimentally the potential of co-registering parameters of cutaneous blood flow and the intrinsic tissue fluorophore fluorescence with the object to determine the stage of lower limb complications of in patients with diabetes mellitus.

The tissue blood perfusion and the tissue fluorescence were assessed in the experimental studies by LDF and FS, respectively (Dremin et al., 2017). Parameters of the blood perfusion and the amplitude of the coenzyme NADH and FAD fluorescence signals were registered simultaneously in the same tissue volume. The single mode laser module with a wavelength of 1,064 nm was used in the laser Doppler channel of the system SPE "LAZMA" Ltd., Russia. For the fluorescence excitation, two radiation sources of 365 and 450 nm were applied. For the assessment of the local regulatory mechanisms of blood flow the local heating stimulation was selected as a provocative action on the blood microcirculation system. A special block "LAZMA-TEST" (SPE "LAZMA" Ltd., Russia) allowed us to establish and maintain the determined local temperature of the investigated area. Specialized software was used to control the unit responsible for thermal functional tests, to observe the real-time course of the experiment and to analyze the recorded parameters.

The experimental study involved 76 patients from the Endocrinology Department of the Orel Regional Clinical Hospital (Oryol, Russia). All patients had type 2 diabetes mellitus. The illness was taking its normal course in 62 patients (Diabetic group 1, mean age 54 ± 10 years). Another 14 people among the patients' group were identified to have a more severe course of the disease (Diabetic group 2, mean age 53 ± 13 years). In each case, the degree of severity was determined based on anamnesis analysis, consultation with the attending physician and presence of trophic disorders: Diabetic group 1 consisted of all patients without trophic ulcers, when Diabetic group 2 included diabetics with the presence of visible trophic ulcers. The control group consisted of 48 healthy volunteers (mean age 46 ± 6 years).

The experimental study was conducted under the established protocol. Each experiment consisted of four stages. The first stage included 4 min registration of microcirculatory blood perfusion under the basal conditions and recording of a pair of fluorescence spectra at the end of the period. This was followed by a consecutive local cold (25◦Ñ) and heat (35 and 42◦Ñ) provocations lasting for 4 min each. The duration of one experiment was approximately 16 min. Time chart of the study protocol is shown in **Figure 11**.

The dorsal surface of the foot was chosen as an investigated area. The probe was secured on a plateau between the 1st and 2nd metatarsal bones (**Figure 12**).

Changes in blood flow of the investigated area of non-glabrous skin are mediated by the vessels themselves and by mechanisms of the sympathetic nervous system (Johnson and Proppe, 2011). The microcirculation assessment in this area may demonstrate better sensitivity to disease severity in comparison with glabrous or non-hairy skin (with arteriovenous anastomoses) where blood perfusion may fluctuate substantially (Thoresen and Walløe, 1980; Vanggaard et al., 2012; Fuchs et al., 2017). All measurements were performed in the supine position. As local pressure of the probe on skin influences strongly on the measurement results of blood perfusion and fluorescence intensity, a particular attention was paid to minimize the

parameter. The volunteers were asked to refrain from food and drink 2 h before the study to exclude the influence of these factors on the change of microhemocirculation and metabolic processes. Room temperature was maintained at a steady 24–25◦C. All volunteers underwent preliminary adaptation to these conditions for at least 10 min. Change of skin temperature was recorded during the study. Since skin temperature of subjects was different in the basal conditions, the studied area was pre-cooled to 25◦C to unify measurements.

#### Data Analysis and Diagnostic Criteria

Averaged amplitudes of fluorescence were recorded at the basal conditions. Normalized intensities of backscattered excitation radiation (AF<sup>460</sup> at 460 ± 10 nm with excitation at 365 nm and AF<sup>525</sup> at 525 ± 10 nm with excitation at 450 nm) were analyzed. The wavelengths were chosen to maximize the probability of detecting NADH and FAD fluorescence signals (Dunaev et al., 2015). The microcirculatory blood perfusion of the investigated tissue was recorded and averaged at each stage of the experiment. Non-parametric methods (Mann–Whitney U-test and Kruskal–Wallis ANOVA test) were used to confirm the reliability of differences in the results. The difference was considered significant when p < 0.01.

Taking into account the division of experimental data into three classes (two diabetic groups and one control), a linear discriminant analysis was used to determine the discriminant functions, which allow stating the desired decision-making rules. Thus, the obtained classifiers allow the newly appearing object to be assigned to one of the above classes by the measuring of its parameters. The quality of the discriminant analysis was assessed using the ROC-curve.

Experimental studies have shown that patients with diabetes mellitus have increased values of normalized fluorescence amplitudes at basal conditions and reduced perfusion response to local heating (up to 35 and 42◦C). It was also noticed that these parameters in the Diabetic group 2 (patients with diabetic trophic ulcers) had significantly difference from the Diabetic group 1 and control group (**Figure 13**).

Statistically significant differences in parameters may be due to a number of reasons. Diabetes can cause tissue hypoxia, which leads to a violation of the aerobic respiration pathway of cells (Ditzel and Standl, 1975; Thangarajah et al., 2009; Xi et al., 2016). In this case, the mitochondrial oxidation of NADH slows down, while the glycolysis pathway of NADH formation activates. Therefore, NADH accumulation may be considered as a sign of tissue hypoxia, and its contribution to the total fluorescence signal as a marker of general oxygen deficiency in tissues.

As was mentioned earlier, AGEs are involved in the mechanisms of diabetes complications development. The rapid formation of intracellular AGEs contributes to the violation of protein function and can serve as an objective marker of glycation in tissues (Gkogkolou and Böhm, 2012). Observed in the study level of skin fluorescence is related to the determined in vitro standard glycation marker HbA1c. The conventional measure of glycation using HbA1c characterizes the glycation

processes that occurs in a short period (around 3 months). While long-term processes in the diabetic skin are reflected by changes in AGE content. In view of the stability of AGEs and the long molecular lifetime of collagen, it might be supposed that the skin fluorescence would be applied as a measure of the total impact of hyperglycaemia in the course of life.

The basal level of perfusion does not have a statistically significant differences between the compared groups (**Figure 13B**). However, the proposed functional thermal tests allow identifying the microvascular disorders. Normally, the increase in perfusion during the local heating occurs mainly due to two mechanisms. As reported, skin heating to 34–35◦C leads to activation of peptidergic nerve fibers by activating heat-sensitive receptors (Stephens et al., 2001). Skin heating to 42◦C leads to vasodilation by releasing NO from the vascular endothelium (Minson et al., 2001; Kellogg et al., 2009). In the case of diabetes, all aspects of the microcirculatory and tissue systems, including vascular endothelial and perivascular nerve fibers, are subject to dysfunction. Therefore, the thermal test is pathogenically justified for the diagnosis of these disorders in diabetes. The reduction of perfusion growth in patients with diabetes mellitus when heated to 35◦C is an objective criterion for impaired function of sensory nerve fibers, which is a component of diabetic neuropathy. The reduction of perfusion growth when heated to 42◦C reflects a deficiency of endothelium-dependent vasodilation mechanisms. Differences in perfusion reaction to various stimuli can be explained by the fact that NO synthesis can be suppressed by the accumulation of AGEs in endothelial cells (Bucala et al., 1991; Stitt et al., 1997; Chakravarthy et al., 1998).

For the synthesis of the decision rule, the analyzed parameters of normalized amplitudes of skin fluorescence and perfusion are proposed. The selected parameters satisfy the principles of statistical independence and the significance of the difference between their values in target groups. Discriminant function was synthesized to provide the best combination of sensitivity and specificity. The lowest level of error was obtained by combining the fluorescence intensity at the excitation wavelength of 365 nm and the average level of skin blood perfusion

at 42◦C. Sensitivity and specificity of the first classification rule (Control group against Diabetic group 1) were 0.92 and 0.90, respectively. Sensitivity and specificity of 0.86 and 0.85, respectively, were obtained for the second classification rule (Diabetic group 1 vs. Diabetic group 2). **Figure 14A** shows the scatter plot of experimental data with applied discriminant lines dividing the experimental points into three groups. Thus, the resulting diagnostic criterion, using two discriminant functions, D1 and D2, and allowing classifying the measured subject to one of three groups, can be written as:

$$\begin{aligned} D1 &= -2.55 - 0.45 \cdot AF\_{460} + 0.23 \cdot I\_{\text{m}}^{42 \text{ } ^\circ \text{C}}, \\ D2 &= -1.80 + 1.16 \cdot AF\_{460} - 0.18 \cdot I\_{\text{m}}^{42 \text{ } ^\circ \text{C}} \end{aligned} \tag{6}$$

**Figure 14A** shows that the shift of the values I<sup>m</sup> <sup>42</sup> ( ◦Ñ and AF460 towards the top-left corner of the diagram characteriszes the deterioration of the patient's condition and increasing the risk of foot ulcers development.

( ◦Ñ and AF460 obtained by LDA method with the applied discriminant functions (6) shown by solid (D2) and dash (D1) lines. The healthy group is shown by red squares, the diabetic group

The ROC-curves calculated for the obtained discriminant functions are shown in **Figure 14B**. The integral characteristic – area under curve (AUC) was used to evaluate and compare the quality of different classification rules. AUC = 0.93 for both found classification rules that means a high level of the classifiers quality.

Thus, the values of skin fluorescence and level of blood perfusion during the local heating test can be used effectively for the analysis the various stages of diabetes, from the initial stage of the disease, to the development of trophic ulcers.

#### Experimental Results and Discussion

The proposed method, based on the use of optical non-invasive diagnostic methods in combination with local heat test, shows the high sensitivity and can be an effective marker of metabolic changes in diabetes mellitus. Both approaches (laser Doppler flowmetry and fluorescence spectroscopy) individually and in combination have excellent diagnostic potential as they can help identify patients at risk of future problems in their lower limbs. The obtained experimental results confirm the prospects of the combined use of these optical non-invasive methods in the detection of biological tissue disorders in type 2 diabetes. In addition, the method has a good potential in the field of evaluation of therapeutic interventions aimed at preventing and reducing the progression of diabetic complications, as well as to minimize their negative consequences.

Expanding the number of clinical trials and improving the methodological support of optical non-invasive diagnostics will allow introducing these methods into clinical practice of the attending physician. The use of this methods of diagnosis of the diabetes complications, which is easily implemented, non-invasive convenient and safe for the patient, will make a significant contribution to the fight against diabetes.

One promising area for further development in this area is the differentiation and evaluation of the contributions of the AGEs, NADH, and FAD to the resulting fluorescence signal. This would allow carrying out studies to identify possible pathways of biological tissue disturbance.

## CONCLUSION

The use of optical non-invasive diagnostic methods has a great potential for the detection of concomitant microcirculation disorders in patients with rheumatic diseases and diabetes. In this review, it was shown that the use of laser Doppler flowmetry, optical tissue oximetry and fluorescence spectroscopy together or separately may have important diagnostic value for the detection of violations, assessment of their severity, as well as for the analysis of the effectiveness of the therapy. The joint application of the considered technologies with the methods of machine learning (discriminant analysis) can be used as additional diagnostic criteria in the field of rheumatology and endocrinology.

Thus, the introduction of optical methods for the assessment of peripheral hemodynamics and metabolic processes in rheumatic diseases and diabetes makes it possible to early preclinical diagnosis of microcirculatory and metabolic disorders, contributing to the improvement of the effectiveness and quality of treatment.

The perspective direction of the development of the optical non-invasive diagnostics method is their realization in the form of the wearable compact wireless devices suitable for long-term monitoring of microcirculation and metabolism parameters accomplished with the development of methodology of their application. The use of compact wearable devices will allow diagnosis of arterio-venular anastomoses and investigation of the compensatory mechanisms and synchronization of blood flow oscillations during functional tests.

# ETHICS STATEMENT

These studies were carried out in accordance with the Declaration of Helsinki principles and approved by the Ethics Committee of the Orel State University. All patients participated in the trials have given the full consent on measurements in written and been informed of their right to discontinue participation at any time.

# AUTHOR CONTRIBUTIONS

AD, SS, and ER designed the study. AZ, VD, IM, and EZ collected and analyzed researched data and drafted the manuscript. AG contributed to the data analysis and preparation of a final version of the article. All authors contributed to the discussion and reviewed and edited manusript.

### FUNDING

SS and ER were partially supported by Russian Science Foundation (Grant No. 18-15-00172). EZ kindly acknowledges

for personal support from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No. 703145 and Academy of Finland, grant No. 318281. AZ (regarding the experimental studies of the blood perfusion and skin temperature dynamics for functional diagnostics of intradermal finger vessels) was funded by the Russian Science Foundation (Grant No. 18-79-00237).

#### REFERENCES


#### ACKNOWLEDGMENTS

We would like to thank all of our volunteers and patients for their contribution to this research project. Special thanks are extended to doctors Khakhicheva L. S., Alimicheva E. A., Masalygina G. I., and Muradyan V. F. of the Orel Regional Clinical Hospital for providing useful advice and help. AG acknowledges the Scottish Informatics & Computer Science Alliance (SICSA) for partly support.



with atherosclerosis: relationship to tissue AGEs. Mol. Med. 3, 617–627. doi: 10.1007/BF03401819


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Zherebtsova, Dremin, Makovik, Zherebtsov, Dunaev, Goltsov, Sokolovski and Rafailov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Oligonucleotide Microarrays Identified Potential Regulatory Genes Related to Early Outward Arterial Remodeling Induced by Tissue Plasminogen Activator

Olga Plekhanova1,2 \*, Yelena Parfyonova1,2, Irina Beloglazova1,2, Bradford C. Berk<sup>3</sup> and Vsevolod Tkachuk1,2

<sup>1</sup> Faculty of Medicine, Lomonosov Moscow State University, Moscow, Russia, <sup>2</sup> National Medical Research Center of Cardiology, Moscow, Russia, <sup>3</sup> Aab Cardiovascular Research Institute, University of Rochester Medical Center, Rochester, NY, United States

#### Edited by:

Alexey Goltsov, Abertay University, United Kingdom

#### Reviewed by:

Alex Bobik, Baker Heart and Diabetes Institute, Australia Dmitry Traktuev, University of Florida, United States

> \*Correspondence: Olga Plekhanova plekhanova@mail.ru

#### Specialty section:

This article was submitted to Vascular Physiology, a section of the journal Frontiers in Physiology

Received: 29 November 2018 Accepted: 08 April 2019 Published: 30 April 2019

#### Citation:

Plekhanova O, Parfyonova Y, Beloglazova I, Berk BC and Tkachuk V (2019) Oligonucleotide Microarrays Identified Potential Regulatory Genes Related to Early Outward Arterial Remodeling Induced by Tissue Plasminogen Activator. Front. Physiol. 10:493. doi: 10.3389/fphys.2019.00493 Constrictive vascular remodeling limiting blood flow, as well as compensatory outward remodeling, has been observed in many cardiovascular diseases; however, the underlying mechanisms regulating the remodeling response of the vessels remain unclear. Plasminogen activators (PA) are involved in many of the processes of vascular remodeling. We have shown previously that increased levels of tissue-type PA (tPA) contributes to outward vascular remodeling. To elucidate the mechanisms involved in the induction of outward remodeling we characterized changes in the expression profiles of 8799 genes in injured rat carotid arteries 1 and 4 days after recombinant tPA treatment compared to vehicle. Periadventitial tPA significantly increased lumen size and vessel area, encompassed by the external elastic lamina, at both one and 4 days after treatment. Among 41 differentially expressed known genes 1 day after tPA application, five genes were involved in gene transcription, five genes were related to the regulation of vascular tone [for example, thromboxane A2 receptor (D32080) or non-selectivetype endothelin receptor (S65355)], and eight genes were identified as participating in vascular innervation [for example, calpain (D14478) or neural cell adhesion molecule L1 (X59149)]. Four days after injury in tPA-treated arteries, four genes, regulating vascular tone, were differentially expressed. Thus, tPA promotes outward arterial remodeling after injury, at least in part, by regulating expression of genes in the vessel wall related to function of the nervous system and vascular tone.

Keywords: tissue type plasminogen activator, outward vascular remodeling, microarray, vascular tone, arterial injury

# INTRODUCTION

Vascular remodeling is one of the most important mechanisms responsible for lumen narrowing in many cardiovascular pathologies (Paneni et al., 2017). Inward arterial remodeling is associated with high cardiovascular mortality, while outward vascular remodeling is considered to be a positive compensatory mechanism that provides adequate blood flow (van Varik et al., 2012). We have

**103**

demonstrated previously that tissue-type plasminogen activator (tPA) attenuated inward arterial remodeling (Parfyonova et al., 2004).

Despite the active interest of many research groups to the problem of regulation of arterial remodeling, highly effective mechanisms are still to be found (Goel et al., 2012). The mechanisms, underlying tPA effects in the vessel wall, stay unclear.

Vascular remodeling involve many processes including cell proliferation, migration, extracellular matrix remodeling, and changes in vascular tone (Alexander and Owens, 2012; van Varik et al., 2012). There are contradictory data about the influence of tPA; it contributes to cell proliferation, but also in the previous study, tPA reduced neointima formation (Parfyonova et al., 2004).

To clarify mechanisms that may lead to the positive outward remodeling induced by tPA, we hypothesized that tPA might affect the expression of genes, which might regulate the remodeling response of the injured artery. To assess the role of gene expression changes in early processes leading to the outward vascular remodeling induced by tPA, we investigated alterations in gene expression profiles of rat common carotids after injury and local treatment with the recombinant tPA.

# MATERIALS AND METHODS

#### Ethics Statement

Animal studies were conducted in accordance with the principles of the Basel Declaration, approved guidelines of the Institutional Animal Care and Use Committee of the National Medical Research Center of Cardiology (Permit No. 385.06.2009) and the Guide for the Care and Use of Laboratory Animals (8th Edition, 2011, The National Academies Press, United States).

### Animal, Experimental Design, Surgery and Tissue Collection

Male Wistar-Kyoto rats (4–5 months old, weight 300–350 g) were obtained from a standard colony. Ketamine anesthesia was used (100 mg/kg body wt) and the balloon catheter (Fogarty 2F) injury of rats left common carotid arteries was performed. tPA (20 nmol/kg, Alteplase, Boehringer Ingelheim Pharma) or 500 µl of vehicle (saline) in Pluronic gel F-127 (BASF) was applied around injured carotids (Plekhanova et al., 2008). At days 1, 4, or 14 after injury the perfusion-fixation was done under anesthesia, the absence of endothelium in the injured area was verified by Evans blue dye (Gangadharan et al., 2001) and paraffin embedding of common carotids was performed (Plekhanova et al., 2000). For microarrays, 1 and 4 days after the operation carotids were rapidly removed and frozen in liquid nitrogen. Total RNA was extracted from left carotid arteries of injured, sham-operated and intact rats (nine animals per group). Microarray analysis was performed (Plekhanova et al., 2008).

#### Morphometry

Toluidine blue staining of cross-sections (10 µm) with subsequent analysis using a Zeiss microscope coupled to a ProgRess-3008 camera (Kontron Elektronik) and PC with an UTHSCSA ImageTool, version 3.0 (San Antonio, TX, United States) (Plekhanova et al., 2000) were performed.

#### Microarray Analysis

Analysis was performed as previously described (Korshunov et al., 2006; Plekhanova et al., 2008). Total RNA was isolated with the use of Qiagen RNAeasy Micro kit, the quality of RNA was confirmed. For each probe vessels from three rats were pooled; three hybridizations were done for each treatment group. Microarray analysis was performed at the University of Rochester Microarray Core Facility in accordance to Affymetrix recommendations using Affymetrix rat genome U34A oligonucleotide microarrays with 8799 known genes and expressed sequence tags (ESTs) (Santa Clara, CA, United States). Streptavidin phycoerythrin stain (SAPE, Molecular Probes) was used. Target hybridization was detected and quantified with a Gene Array Scanner (Hewlett Packard/Affymetrix). "Array performance" was assessed. Hybridization efficiency and sensitivity was controlled using control transcripts spiked into the hybridization cocktail. Inter-array variability was assessed by Microarray Analysis Suite 5.0. For data normalization global scaling with target intensity of 500 across all probe sets was used (Brooks et al., 2002). To lower the rates of overall errors and the estimated false detection the invariant set normalization method and the model-based expression indexes for perfect match (PM)-only arrays using DNA-Chip Analyzer (dChip) software package were applied (Li and Hung Wong, 2001). The "Present" percentage in all arrays was >40%. The percentage of "Array-outlier" in all arrays was not more than 5%, showing the absence of serious contamination or hybridization problems. To exclude "Absent" genes filtering was performed limiting the analyses to targets, which were "present" in more than 20% of arrays. A minimal inter-array variability (SD ≤ 2%) of the mean signal (with 5% of signals trimmed from both the high and low ends) was observed across the 4458 targets after the presence/absence filtering. Then to identify statistically significant changes in gene expression the significance analysis of microarrays (SAM) software was applied (Welle et al., 2002). Comparison with parametric ANOVA and Welch's approximate t-test for non-equal variances were also performed (GeneSpring Software). To get very robust results we accepted as significant only genes verified by both analyses (Hoffmann et al., 2002).

#### Software

Microarray Suite 5.0 was used with the default parameters. DNA-Chip Analyzer software was kindly provided by Drs. Cheng Li and Wing Hung Wong, Computational Biology Lab, Department of Biostatistics, Harvard School of Public Health. SAM was kindly provided by the Department of Statistics, Stanford University. GeneSpring 5.0 Software was from Silicon Genetics.

## Quantitative Polymerase Chain Reaction (qPCR)

Quantitative RT-PCR analyses were performed using ABI Prism 7900HT sequence detection system (Applied Biosystems).

Double-stranded cDNA template preparation and purification were performed with Ambion MessageAmp aRNA kit. The qPCR primers and Master Mix from RT2 Real-Time Gene Expression Assay kits (SuperArray) were used. The effective hybridization was testified by the expression of 7 housekeeping genes including ribosomal proteins L18 (M20156), L18a (X14181) and S9 (X66370), alpha-tubulin (V01227), thymosin beta-4 (M34043), ubiquitin (D17296), beta-actin (V01217), which consistently gave positive signals through all samples analyzed.

#### Statistical Analysis

Results are represented as mean ± SEM. Jandel SigmaStat (t-test or ANOVA for multiple comparisons to limit false discovery rate, as appropriate) was used for statistical analysis. P < 0.05 was considered as statistically significant.

# RESULTS

## tPA and Vessel Remodeling

To assess tPA influence on the injured artery remodeling, we measured the area of lumen, the area, encompassed by the external elastic lamina (EEL), and intima-plus-media area [termed intima-media thickness (IMT)]; main parameters reflecting vessel remodeling. Periadventitial application of tPA increased EEL area (P < 0.05; **Figure 1**) at 1 day, then at 4 days there were significant increases in EEL and lumen areas, and a decrease in IMT (P < 0.05; **Figure 1**). Fourteen days after tPA treatment, there was still a decrease in IMT, but the other areas were the same (P < 0.05; **Figure 1**). These effects contrasted markedly with the effects of uPA that induced lumen narrowing, decreases in EEL and increased IMT after injury (Plekhanova et al., 2001, 2008). These data show that the arterial remodeling response stimulated by tPA induces lumen enlargement (**Figure 1**), increased EEL and decreased IMT, although the decrease in IMT only lasted 14 days.

### Differential Gene Expression Profiling Induced by tPA

To clarify the underlying mechanisms ensuring the vascular remodeling effects of tPA, the in vivo microarray analysis of gene expression was carried out. The expression of 8799 genes was investigated 1

and 4 days after the combination of balloon arterial injury plus treatment with recombinant tPA or control (gel plus saline).

One day after the operation in arteries treated with tPA the expressions of 95 transcripts (41 known genes, **Supplementary Table 1A**) were significantly different to those in control samples (P ≤ 0.05; confirmed by two statistical tests). Wherein amongst known genes 23 (56.1%) were downregulated, when 18 (43.9%) were up-regulated after tPA treatment compared to control. At day 4, 53 mRNAs (28 known genes, **Supplementary Table 1B**) were significantly different in tPA group compared to gel-saline group. Among these known genes 5 (17.9%) were down-regulated and 23 (82.1%) were up-regulated compared to control.

In the group of genes differentially expressed after tPA treatment, a few groups of genes participating in the main processes of vascular remodeling including cellular migration, apoptosis and proliferation and inflammation were identified (**Supplementary Tables 1A,B**).

#### tPA and Differential Expression of Nervous System – Related Genes

Among genes differentially expressed 1 day after injury and tPA treatment, eight genes were related to nervous system function (**Table 1**). Among these genes, seven were up-regulated, and one was down-regulated. At day four after injury the expression of one gene with the intended relevance to nervous system was increased in arteries treated with tPA. Changes in genes expressions of ≥1.5 fold were observed (**Supplementary Figure 1**). Thus, microarray gene analysis showed that tPA may contribute to the outward vascular remodeling at least partially through local changes in the expression of nervous system genes after arterial injury.

## tPA – Induced Differential Regulation of Vascular Tone Related Genes

Genes involved in regulation of vascular tone were discovered to be differentially expressed in arteries both one and 4 days after injury and tPA application (**Table 1**). One day after injury and tPA treatment, one gene was up-regulated and one gene was down-regulated; while at 4 days after injury, three genes were up-regulated and 1 gene was down-regulated. Changes in genes expressions of 1.5–3.0 folds were observed (**Supplementary Figure 2**). These findings suggest the importance of tPA-mediated changes in vascular tone in vascular remodeling.

#### DISCUSSION

We demonstrated three major results regarding the role of tPA in vascular remodeling after injury. First, exogenous tPA significantly enlarged arterial lumen area at 1 as well as 4 days

TABLE 1 | Transcription profiles of genes related to nervous system and vascular tone regulation differentially expressed in rat carotid arteries treated with tPA compared to control vessels, which only received Pluronic gel/saline, 1 day and 4 days after balloon injury.


after balloon injury. Second, we discovered increases in EEL area at both one and 4 days. Finally, tPA decreased IMT at 4 and 14 days. These results show that tPA has a significant effect on vessel structure that appears to be beneficial. Our data are consistent with previous studies where local delivery of tPA via Dispatch catheters, followed by continuous intravenous infusion of tPA for 3 days, prevented intimal hyperplasia after angioplasty (Kanamasa et al., 2001). In contrast, an adenoviral construct expressing tPA enhanced neointima formation after angioplasty in a rabbit model (Hilfiker et al., 2001). Another study in cholesterol-fed rabbits showed that subcutaneous recombinant tPA inhibited the effect of balloon injury on luminal area compared with controls, but this effect did not reach statistical significance; the potential problem with the delivery method was discussed (Alexander et al., 2004). Taken together, these results indicate the importance of the model and the method used for tPA administration to evaluate its effects on arterial remodeling.

The gene array data shown here provide information regarding two likely mechanisms of action of tPA in the vessel wall: nerve function and vascular tone. Changes in nervous system–related gene expression after tPA treatment of injured arteries suggest that tPA-induced increases of vascular lumen size can be mediated at least in part through its effects on innervation that is consistent with the previous findings (Stenmark et al., 2011; Chistiakov et al., 2015). Some of these genes may affect vasodilatation. For example, NMDA receptor through nitric oxide (Faraci and Breese, 1993; Hama-Tomioka et al., 2012) or substance P receptor may elicit a neurogenic vasodilatation (Galligan et al., 1990; Charkoudian, 2010). Future work will be necessary to define precisely the roles of neurotransmitters in arterial remodeling. Of interest among the other nerve function genes, calpain was found to be highly expressed in the vasculature and is involved in many processes important for vascular remodeling (Letavernier et al., 2008). The involvement of its signaling in vascular remodeling regulation is extensively investigated nowadays (Kovacs and Su, 2017). Also smooth muscle cells neurotrophin-3 was shown to be implicated in vessel injury acting in an autocrine manner through its receptor TrkC (Donovan et al., 1995).

The microarray study also revealed in tPA-treated arteries changes in the expression of four genes related to vascular tone regulation, such as endothelin receptor or angiotensin receptor (Watanabe et al., 2005; Rapoport and Merkus, 2017). One of those genes - corticotrophin releasing factor receptor (U53486) is also involved in nervous system function. Considering that perivascular nerves may also affect vascular tone and induce vasodilatation (as well as tPA itself possesses vasoactive properties) (Heyman et al., 2004; Nassar et al., 2004), we may hypothesize that increases in lumen and EEL areas early after injury and tPA application may reflect changes in vascular tone, though it was not pronounced at the later stages due to the nature of vascular tone response.

The most novel finding here is the discovery that changes in genes related to nervous system function as well as vascular tone regulation, were highly significant at early time points, suggesting a possible connection between these two pathways leading to early outward vessel remodeling. There are many examples of relationships between nerves and vessels, especially regarding vascular tone, such as the release of vasoactive intestinal peptide by nerves to promote vasodilation in the gut.

In summary, the present data further support an important role for tPA as a negative regulator of IMT in response to balloon injury, as well as a mediator of favorable arterial remodeling. Furthermore, our genetic analysis reveal a novel role for the interaction between genes involved in the nervous system and those in vascular tone regulation suggesting new approach to limit IMT after injury.

# ETHICS STATEMENT

Animal studies were conducted in accordance with the principles of the Basel Declaration, approved guidelines of the Institutional Animal Care and Use Committee of the National Medical Research Center of Cardiology (Permit No. 385.06.2009) and the Guide for the Care and Use of Laboratory Animals (8th Edition, 2011, The National Academies Press, United States).

# AUTHOR CONTRIBUTIONS

OP carried out the experiments and wrote the manuscript. YP and IB provided consultancy. BB supervised the project and wrote the manuscript. VT supervised the project.

# FUNDING

This study was supported by National Institutes of Health grant HL-62836 to BB, by NHLBI funds for exchanges of scientists under the auspices of the United States–Russia Joint Agreement in Cardiopulmonary Research, by NATO LST. CLG.979209, by Grants of the Russian Science Foundation (No. 14-24-00086 in part of gene expression profiling induced by tPA; No. 19-75-30007 in part of vessel and tissue remodeling processes).

# ACKNOWLEDGMENTS

The authors are grateful for assistance of Drs. Pavel Bashtrykov and Oleg Viborov with animal procedures and Marina Solomatina with morphometry.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019. 00493/full#supplementary-material

# REFERENCES

fphys-10-00493 April 27, 2019 Time: 15:33 # 6


increased carotid intima–media thickening. Arterioscler. Thromb. Vasc. Biol. 26, 295–300. doi: 10.1161/01.ATV.0000196544.73761.82


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Plekhanova, Parfyonova, Beloglazova, Berk and Tkachuk. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

digital media

of impactful research

article's readership