A Novel MRA-Based Framework for Segmenting the Cerebrovascular System and Correlating Cerebral Vascular Changes to Mean Arterial Pressure

Featured Application: This non-invasive methodology could be used to detect alterations in the cerebrovasculature by analyzing MRA images, which would assist clinicians to optimize medical treatment plans of HBP. Abstract: Blood pressure (BP) changes with age are widespread, and systemic high blood pressure (HBP) is a serious factor in developing strokes and cognitive impairment. A non-invasive methodology to detect changes in human brain’s vasculature using Magnetic Resonance Angiography (MRA) data and correlation of cerebrovascular changes to mean arterial pressure (MAP) is presented. MRA data and systemic blood pressure measurements were gathered from patients (n = 15, M = 8, F = 7, Age = 49.2 ± 7.3 years) over 700 days (an initial visit and then a follow-up period of 2 years with a ﬁnal visit.). A novel segmentation algorithm was developed to delineate brain blood vessels from surrounding tissue. Vascular probability distribution function (PDF) was calculated from segmentation data to correlate the temporal changes in cerebral vasculature to MAP calculated from systemic BP measurements. A 3D reconstruction of the cerebral vasculature was performed using a growing tree model. Segmentation results recorded 99.9% speciﬁcity and 99.7% sensitivity in identifying and delineating the brain’s vascular tree. The PDFs had a statistically signiﬁcant correlation to MAP changes below the circle of Willis ( p -value = 0.0007). This non-invasive methodology could be used to detect alterations in the cerebrovascular system by analyzing MRA images, which would assist clinicians in optimizing medical treatment plans of HBP.


Introduction
High blood pressure (HBP) affects approximately 1 in 3 adults in the USA. HBP is a primary or a contributing cause of mortality in about 410,000 adults each year with associated healthcare costs of $46 billion [1]. High sodium intake [2], chronic stress [3], and renal dysfunction [4] are the primary causes for HBP. Clinical studies suggest that hypertension development is correlated to changes in the vascular structure of human brains [5]. These cerebrovascular changes are hypothesized to be a significant contributor to strokes, brain lesions, cerebral ischemic injury, dementia and cognitive impairment [5][6][7].
Currently, HBP is diagnosed and medically managed when systemic BP measurements using sphygmomanometer are greater than 140/90 mmHg. However, BP measurement via sphygmomanometer cannot quantify cerebrovascular structural changes that can increase the risk of cerebral adverse events. Recently, some studies hypothesized that vascular structural changes in human brains may occur prior to the elevation of systemic BP rather than cerebrovascular damage due to sustained exposure to HBP [8][9][10][11]. Therefore, quantification of cerebral vascular alterations could help identify and stratify patients who have potential of cerebral adverse events, potentially enable managing appropriate medical plans prior to the onset of systemic hypertension in conjunction with other cognitive tests and optimize medical management of HBP patients.
Magnetic Resonance Imaging (MRI) technique has been traditionally used in the quantification of organ structural changes. In the literature, MRI scanning has been used for volumetric measurement of the ventricular cavities and myocardium [12] and to determine intravascular pressures from magnetic resonance (MR) velocity data in large vessels like the aorta or pulmonary artery [13]. MRA scanning has been used to quantify measurements of the flow in the collateral arteries of patients that have occlusions in internal carotid artery [14]. To the best of our knowledge, neither MRA nor MRI has been utilized for the estimation of vascular pressure changes in the brain. Detection of cerebrovascular changes in MRA images has not been accomplished due to the lack of accurate segmentation algorithms that can extract tiny blood vessels in human brains (in comparison to aorta or pulmonary arteries) from the surrounding soft tissue. Further, there are no methods to quantify cerebrovascular structural changes or to correlate them to changes in mean arterial pressure, MAP from MRI/MRA imaging.
The goal of this manuscript is to develop a new framework that detects the potential changes in cerebrovascular structure by first introducing a novel automatic segmentation algorithm that delineates the cerebrovascular system from MRA data, and then estimate the change in cerebral vascular diameters to demonstrate proof-of-concept of correlation between cerebrovascular structural changes to MAP.

Materials and Methods
In this section, patient demographics and details about the proposed methodology and data analysis are presented.

Patient Demographics
This work has been approved by the Institutional Review board (IRB) at Pittsburgh University. MRA scans and corresponding blood pressure measurements were acquired from subjects (n = 15, M = 8, F = 7, age = 49.2 ± 7.3) during a 700 day study and were analyzed retrospectively. Participants in this study were carefully selected to represent a wide range of BP changes over the 700 day period and the MRA data were analyzed blinded to the corresponding BP measurements. Participants were chosen in accordance to some exclusion criteria: (1) The use of prescribed medication for hypertension and psychotropic; (2) medical conditions including pregnancy, ischemic coronary artery disease, chronic kidney disease (creatinine > 1.2 mg/dL), chronic liver disease, diabetes mellitus (fasting blood glucose > 125 mg/dL), or cancer (treatment < 12 months); (3) neuropsychiatric conditions including serious head injury, multiple sclerosis, epilepsy, stroke, major mental illness, and brain tumor. A 3T Trio TIM whole-body (Siemens Medical Solutions, Erlangen, Germany) scanner with a 12-channel phased-array head coil was used to acquire the MRA imaging data. Blood pressure and MAP values were calculated from an average of four sphygmomanometer readings taken during two visits before the MRA scanning [9]. First, participants were asked to sit down with their back and arm supported for at least 5 min, and then trained assistants measured BP twice, with 60 s separation time. Assistants used the auscultatory technique with cuff size based on arm circumference to measure BP. An average of the four readings across the two visits was calculated as the resting BP. Each MRA scan consists of 3D multi-slab high resolution images with 160 slices, thickness of 0.5 mm, resolution of 384 × 448. Moreover, the voxel size is 0.6 × 0.6 × 1.0 mm, the flip angle is 15 degrees, the repetition time is 21 ms, and the echo time is 3.8 ms. The scan takes 7 min 39 s and did not require any contrast.
The subjects had an average day 0 systolic pressure of 122 ± 6.9 mmHg, an average day 0 diastolic pressure of 82 ± 3.8 mmHg, at an average day 700 systolic pressure of 118.9 ± 12.4 mmHg, and an average day 700 diastolic pressure of 79.9 ± 11.0 mmHg (i.e., mean systolic pressure remained comparable over time though some individuals increased in pressure and some decreased or stayed the same).

Data Analysis
The analysis of patient MRA data consists of five key steps ( Figure 1): (1) manual segmentation of training slices to identify ground truth, (2) automatic segmentation for all slices to delineate the blood vessels from the surrounding soft tissue by combining the segmented ground truths with the Linear Combination of Discrete Gaussians (LCDG) models for grey level distribution, (3) voxel matching for obtaining temporal subtraction images to enhance the ability to see cerebrovascular change via a distance map created to quantify the change in patients between day 0 and day 700, (4) generation of a probability distribution function (PDF) which describes the distribution of pixel distances from vascular edges and is used to statistically correlate to BP and (5) estimation of cumulative distribution function (CDF) to observe the summated probability of cerebrovascular changes in the same patient from day 0 and day 700. Each MRA scan consists of 3D multi-slab high resolution images with 160 slices, thickness of 0.5 mm, resolution of 384 × 448. Moreover, the voxel size is 0.6 × 0.6 × 1.0 mm, the flip angle is 15 degrees, the repetition time is 21 ms, and the echo time is 3.8 ms. The scan takes 7 min 39 s and did not require any contrast. The subjects had an average day 0 systolic pressure of 122 ± 6.9 mmHg, an average day 0 diastolic pressure of 82 ± 3.8 mmHg, at an average day 700 systolic pressure of 118.9 ± 12.4 mmHg, and an average day 700 diastolic pressure of 79.9 ± 11.0 mmHg (i.e., mean systolic pressure remained comparable over time though some individuals increased in pressure and some decreased or stayed the same).

Data Analysis
The analysis of patient MRA data consists of five key steps ( Figure 1): (1) manual segmentation of training slices to identify ground truth, (2) automatic segmentation for all slices to delineate the blood vessels from the surrounding soft tissue by combining the segmented ground truths with the Linear Combination of Discrete Gaussians (LCDG) models for grey level distribution, (3) voxel matching for obtaining temporal subtraction images to enhance the ability to see cerebrovascular change via a distance map created to quantify the change in patients between day 0 and day 700, (4) generation of a probability distribution function (PDF) which describes the distribution of pixel distances from vascular edges and is used to statistically correlate to BP and (5) estimation of cumulative distribution function (CDF) to observe the summated probability of cerebrovascular changes in the same patient from day 0 and day 700. Manual Segmentation of Training Slices: MRA data from a patient consists of 160 MRA slices. Every tenth slice was manually segmented to extract the blood vessels from surrounding tissue using Adobe Photoshop (Adobe Systems, San Jose, CA, USA). This methodology allows for delineation of the in-plane blood vessel from the surrounding tissue at a pixel level accuracy where the largest limitation is the resolution of the MRI machine itself. The manually segmented training binary (black for surrounding tissue and white for target vasculature) slices are referred to as ground truths (GT) as the images are correct and free from artifacts or noise ( Figure 2). The manual segmentation of select slices was used for the initialization and optimization of the segmentation algorithm, which was subsequently used for segmenting all obtained slices. Manual Segmentation of Training Slices: MRA data from a patient consists of 160 MRA slices. Every tenth slice was manually segmented to extract the blood vessels from surrounding tissue using Adobe Photoshop (Adobe Systems, San Jose, CA, USA). This methodology allows for delineation of the in-plane blood vessel from the surrounding tissue at a pixel level accuracy where the largest limitation is the resolution of the MRI machine itself. The manually segmented training binary (black for surrounding tissue and white for target vasculature) slices are referred to as ground truths (GT) as the images are correct and free from artifacts or noise ( Figure 2). The manual segmentation of select slices was used for the initialization and optimization of the segmentation algorithm, which was subsequently used for segmenting all obtained slices.
Automatic Segmentation: One of the most challenging issues relating to common computer-assisted diagnostics is the segmentation of accurate 3D cerebrovascular system information from MRA images. Our approach was to rapidly and accurately extract the blood vessel data by defining the probability models for all regions of interest within the statistical approach and not predefining the probability models [15][16][17]. For each MRA slice, the empirical gray level distribution was closely approximated with an LCDG. Then, it was divided into three individual LCDGs, one for every region of interest associated with each of the following dominant modes: gray brain tissues, bright blood vessels, and darker bones and fat. The identified models specify an intensity threshold to extract blood vessels in that slice. A 3D connectivity filter was then applied on the extracted voxels (voxel = volume x element; a representation for a 3-D pixel) to select the desired vascular tree. This method results in higher precision region models with higher segmentation accuracy compared to other methods [16]. Automatic Segmentation: One of the most challenging issues relating to common computer-assisted diagnostics is the segmentation of accurate 3D cerebrovascular system information from MRA images. Our approach was to rapidly and accurately extract the blood vessel data by defining the probability models for all regions of interest within the statistical approach and not predefining the probability models [15][16][17]. For each MRA slice, the empirical gray level distribution was closely approximated with an LCDG. Then, it was divided into three individual LCDGs, one for every region of interest associated with each of the following dominant modes: gray brain tissues, bright blood vessels, and darker bones and fat. The identified models specify an intensity threshold to extract blood vessels in that slice. A 3D connectivity filter was then applied on the extracted voxels (voxel = volume x element; a representation for a 3-D pixel) to select the desired vascular tree. This method results in higher precision region models with higher segmentation accuracy compared to other methods [16].
Adapting the Expectation-Maximization (EM) technique to the LCDG allows for precise identification of the LCDG model which included the number of its components (negative and positive) [18] and for identification of a continuous LCDG-model that contains the probability distribution.
An expected log-likelihood was used as a criterion for model identification in [16]. Consider X = (Xs: s = 1,..., S) to be denoting a 3D MRA slice that contains S co-registered 2D slices Xs = (Xs (i,j): (i,j) ∈ R; Xs (i,j) ₊∈ Q). R and Q = {0, 1,..., Q − 1} represent a rectangular arithmetic lattice that supports the 3D slice and a finite set of Q-ary intensities (gray levels), respectively. Consider Fs = (fs (q): q ∈ Q; ∑q ∈ Q fs (q) = 1), where q is the gray level, to be an empirical marginal probability distribution for gray levels of the MRA slice Xs.
As explained in [18], each slice is considered as a K-modal image with a known number K of the dominant modes related to the regions of interest. For the segmentation of the slice by modes separating, an estimation of the individual probability distributions of the signals associating each mode from Fs is necessary. Fs is closely approximated with LCDG opposing conventional mixture of Gaussians, one per region, or slightly more flexible mixtures involving other simple distributions, one per region. The image LCDG is then divided into sub-models that are related to each dominant mode [19][20][21].
A discrete Gaussian distribution is defined on the set of integers (gray levels) Q = {0,1,...,Q − 1} by the probability mass function: where the parameter θ = (μ, σ), and Φ is the CDF of a normal distribution with mean μ and variance σ 2 . Then the LCDG with Cp positive components and Cn negative components, such that Cp ≥ K, has the probability mass function Adapting the Expectation-Maximization (EM) technique to the LCDG allows for precise identification of the LCDG model which included the number of its components (negative and positive) [18] and for identification of a continuous LCDG-model that contains the probability distribution.
An expected log-likelihood was used as a criterion for model identification in [16]. Consider X = (X s : s = 1, . . ., S) to be denoting a 3D MRA slice that contains S co-registered 2D slices X s = (X s (i,j): (i,j) ∈ R; X s (i,j) + ∈ Q). R and Q = {0, 1, . . ., Q − 1} represent a rectangular arithmetic lattice that supports the 3D slice and a finite set of Q-ary intensities (gray levels), respectively. Consider F s = (f s (q): q ∈ Q; ∑q ∈ Q f s (q) = 1), where q is the gray level, to be an empirical marginal probability distribution for gray levels of the MRA slice X s .
As explained in [18], each slice is considered as a K-modal image with a known number K of the dominant modes related to the regions of interest. For the segmentation of the slice by modes separating, an estimation of the individual probability distributions of the signals associating each mode from F s is necessary. F s is closely approximated with LCDG opposing conventional mixture of Gaussians, one per region, or slightly more flexible mixtures involving other simple distributions, one per region. The image LCDG is then divided into sub-models that are related to each dominant mode [19][20][21].
A discrete Gaussian distribution is defined on the set of integers (gray levels) Q = {0,1, . . ., Q − 1} by the probability mass function: where the parameter θ = (µ, σ), and Φ is the CDF of a normal distribution with mean µ and variance σ 2 . Then the LCDG with C p positive components and C n negative components, such that C p ≥ K, has the probability mass function p w,Θ (q) = ∑ C p r=1 w p,r ψ q θ p,r − ∑ C n l=1 w n,l ψ(q|θ n,l ). (1) The weights w = (w p ,1, . . ., w p ,C p , w n ,1, . . ., w n , C n ) are restricted to be all nonnegative and to satisfy ∑ C p r=1 w p,r − ∑ C n l=1 w n,l = 1.
In general, valid probabilities are nonnegative: p w,Θ (q) ≥ 0 for all q ∈ Q. This implies that the probability distributions only make use of a valid subset of all the LCDGs in (1), which can have negative components p w,Θ (q)< 0 for some q ∈ Q.
Our aim is finding a K-modal probability model which approximates closely the unknown marginal distribution of the gray level. Consider F s , its Bayesian estimate F is as follows [22]: f(q) = (|R|f s (q) + 1)/(|R| + Q), and the intended model should maximize the expected log-likelihood of the statistically independent empirical data with the parameters of the model: The segmentation algorithm basic steps are following [16]: (1) For every slice X s , s = 1,. . . . . ., S, (a) First is to gather the marginal empirical probability distribution F s of gray levels.
Find a starting LCDG model which is nearing F s by using the initialization algorithm to approximate the values of C p −K, C n , and the parameters w, Θ (weights, means, and variances) of the negative and positive discrete Gaussians (DG Use intensity threshold t to extract the voxels of the blood vessels in the MRA slice, which separates their LCDG-sub model from the background. (2) Remove the artifacts from the extracted voxels whole set with a connection filter which chooses the greatest connected tree system built by a 3D growing algorithm [23]. Algorithm 1 summarizes the adopted segmentation approach.

Algorithm 1. Segmentation Approach Main Steps.
For every slice X s , the following steps were completed: 1. LCDG Initialization: • Find the marginal empirical probability distribution of gray levels F s . • Estimate C p − K, C n , W, and Θ of the positive and negative DGs.

•
Find the initial LCDG-model that approximates F s .

LCDG Refinement:
• Fixing C p and C n , refine the LCDG-model with the modified EM algorithm by manipulating other parameters.

Initial Segmentation:
• Divide the final LCDG-model into K sub models by minimizing the expected errors of misclassification.

•
Select the LCDG-sub model that has the largest mean value to be the model of the wanted vasculature.

•
Use the intensity threshold t to extract the voxels of the blood vessels in the MRA slice, separating their LCDG-sub model from the background.

Final Segmentation:
• Remove the artifacts from the extracted voxels whole set with a connection filter which chooses the greatest connected tree system built by a 3D growing algorithm.
This procedure aims to decipher threshold for every MRA image which will enable the complete extraction of the bright blood vessels while removing the darker unwanted tissue and also separating surrounding non-vasculature tissue that may be of similar brightness and along the same boundaries.
Step 1b's initialization creates the LCDG with the non-Appl. Sci. 2021, 11, 4022 6 of 14 negative starting probabilities p w,θ (q).The refinement in 1c increases the likelihood, but the probabilities continue to be non-negative. The experiments presented in [16] show the opposite situations were never met.
Accuracy of the automatic segmentation is evaluated by calculating total error compared to the ground truths. True positive (TP), false positive (FP), true negative (TN), and false negative (FN) segmentations are measured for evaluation.
In Figure 3, if C is the segmented region, G is the ground truth, and R represents the entire image frame, then the TP = |C ∩ G|, the TN = |R − C ∪ G|, the FP = |C − C ∩ G|; and the FN = |G − C ∩ G|. The total error ε is given in [24] as ε = (FN + FP)/(TP + FN) = (FN + FP)/G.

Final Segmentation: •
Remove the artifacts from the extracted voxels whole set with a connection filter which chooses the greatest connected tree system built by a 3D growing algorithm.
This procedure aims to decipher threshold for every MRA image which will enable the complete extraction of the bright blood vessels while removing the darker unwanted tissue and also separating surrounding non-vasculature tissue that may be of similar brightness and along the same boundaries.
Step 1b's initialization creates the LCDG with the non-negative starting probabilities pw,θ (q).The refinement in 1c increases the likelihood, but the probabilities continue to be non-negative. The experiments presented in [16] show the opposite situations were never met.
Accuracy of the automatic segmentation is evaluated by calculating total error compared to the ground truths. True positive (TP), false positive (FP), true negative (TN), and false negative (FN) segmentations are measured for evaluation.
In Figure 3, if C is the segmented region, G is the ground truth, and R represents the entire image frame, then the TP = |C ∩ G|, the TN = |R − C ∪ G|, the FP = |C − C ∩ G|; and the FN = |G − C ∩ G|. The total error ε is given in [24] as ε = (FN + FP)/(TP + FN) = (FN + FP)/G. Voxel Matching: A voxel is an array of volume elements that constitute a notional 3D space. A 3D affine registration is used to handle the pose, orientation, and the data spacing changes and other scanning parameter changes between day 0 and day 700 [25]. In this step, the determined Euclidian radii are converted into diameter values. The output is then converted into a distance map.
Generation of Probability Distribution Function and Validation: The EM-based technique is adapted to the LCDG-model and the distribution of pixel distances is extracted from the distance map to calculate the probability distribution of the cerebrovascular changes. The PDF marks the distribution of white pixels as a true value and black pixels being ignored for the data set. The diameters of the blood vessels are determined by estimating Euclidian center point distances from the edge of a vessel. The data points in the generated PDFs are then extracted and compared to the blood pressure data using statistical analysis.
Calculation of Cumulative Distribution Function: The integral of the PDF is used to generate the CDF as follows: CDF (FX) of a random variable (X) is determined from its Voxel Matching: A voxel is an array of volume elements that constitute a notional 3D space. A 3D affine registration is used to handle the pose, orientation, and the data spacing changes and other scanning parameter changes between day 0 and day 700 [25]. In this step, the determined Euclidian radii are converted into diameter values. The output is then converted into a distance map.
Generation of Probability Distribution Function and Validation: The EM-based technique is adapted to the LCDG-model and the distribution of pixel distances is extracted from the distance map to calculate the probability distribution of the cerebrovascular changes. The PDF marks the distribution of white pixels as a true value and black pixels being ignored for the data set. The diameters of the blood vessels are determined by estimating Euclidian center point distances from the edge of a vessel. The data points in the generated PDFs are then extracted and compared to the blood pressure data using statistical analysis.
Calculation of Cumulative Distribution Function: The integral of the PDF is used to generate the CDF as follows: CDF (F X ) of a random variable (X) is determined from its PDF (f X ) using F X (x) = x −∞ f x (t)dt. CDF shows the summation of the probability that a blood vessel will take a value less than or equal to a vascular diameter value, that defines the blood vessel diameter average in every slice. It shows the cumulative distribution of the PDF with an upper limit of 1. The more quickly the CDF line approaches 1, the more certain that the diameter of the blood vessel is smaller compared to a CDF that takes longer to approach 1. This is illustrated in the results section.

Statistical Analysis
Data were statistically analyzed using R-software (version 3.30) by The R Foundation for Statistical Computing, Vienna, Austria. A mixed effects linear model was used to test the relationship of MRA data with clinical BP measurements. Brain slices were separated into upper compartment (above circle of Willis) and lower compartment (below circle of Willis) to determine correlation with clinical BP readings. Circle of Willis, near the brain base, is where the intracranial cerebral arteries take off from and give rise to progressively smaller vessels [5]. The BP measurements were combined into a single value, the estimated mean arterial pressure (MAP = (2 × DBP + SBP)/3), which was a covariate in the model. Also included in the model were patient's gender, age, and a random intercept. The dependent variable was the mean of Euclidean distance map over the whole vascular tree within each compartment. (Two separate models were fit to the upper and lower compartments.) Statistical significance of fixed effects in the fitted models was determined using likelihood ratio chi-square tests.

3D Reconstruction of the Cerebral Vasculature
A growing tree model that eliminates any unwanted segmented voxels by choosing the greatest connected vascular tree system, coupled with a smoothing algorithm, was used to generate a 3-D model based on segmented slices [23]. An example for the resultant vascular system is visualized and illustrated in the results section. Table 1 shows the specificity and sensitivity of the segmentation algorithm. The automatically segmented slices for all 15 patients were compared to GTs (manually segmented) to calculate accuracy of the algorithm (Figure 4). A cumulative sensitivity of 0.997 ± 0.008 (sensitivity range = 0.969 to 1) and the cumulative specificity of 0.9998 ± 0.0001 (specificity range = 0.9994 to 1) were recorded. The manually segmented training slices (every 10th slice) were excluded in the accuracy calculations of the proposed segmentation approach.  The results of the analysis of the linear mixed effects models (Table 2) revealed that MAP is inversely related to the mean vessel diameter below the circle of Willis (p = 0.0007). The mean diameter of vessels below the circle of Willis was not found to vary significantly with age of the patient or the gender of the patient. Above the circle of Willis, the mean diameters of vessels showed statistically significant decrease with age (p = 0.0005).  with age of the patient or the gender of the patient. Above the circle of Willis, the mean diameters of vessels showed statistically significant decrease with age (p = 0.0005). In the analysis, 13 out of 15 patients showed significant correlation between MAP and the diameters indicated via CDF. Out of the 13 patients that showed CDF correlation with MAP, two example patients (A and B) are shown with the two patients (C and D) where the correlation between CDF and MAP was not found ( Figure 5). Patient C had a shift in CDF in opposition to the MAP change, and patient D had a larger shift in CDF compared to the MAP change ( Figure 5C,D, Table 3). The 3D cerebrovascular model reconstruction of patients C and D indicated significant vascular changes between day 0 and day 700 ( Figure 6).     Figure 5. Sample patient CDFs demonstrating the temporal changes from day 0 to 700. The graphs indicate the probability that blood vessels may be of a certain diameter or less. (A-D) represents 4 different patients from our dateset.

Discussion
The average cumulative segmentation algorithm had a sensitivity of 0.997 ± 0.008 and a specificity of 0.9998 ± 0.001. This high level of accuracy demonstrates the benefit of using a manual input to initialize automatic segmentation. Using manual segmentation alone would be too time intensive to be used in a practical healthcare setting, while utilizing only an automatic segmentation approach would not provide sufficient segmentation accuracy to delineate and quantify the diameters of the smaller arteriolar (<10 micrometers) cerebral blood vessels. The proposed segmentation algorithm combines the accuracy of manual segmentation with the benefit of automated and less time intensive approach and provides segmentation with a high degree of accuracy while also minimizing the required time and effort. Almost all the state-of-the-art algorithms that were used to detect blood vessels were only suitable for healthy blood vessels as they usually assume the vascular linearity and/or circular cross-sections which is not the case in pathological vessels [23]. The proposed segmentation algorithm, however, did not impose any assumptions on blood vessels, and it can delineate both healthy and pathological vessels efficiently. However, the focus in this study was the investigation of the correlation between hypertension development and changes in cerebral vascular diameter. So, the change in a vessel diameter was the only studied diagnostic factor so far. In future work, we would introduce more investigation on more diagnostic factors which in turn would help in diagnosing various vascular-based disease such as strokes.
The high degree of sensitivity and specificity of our approach in accurately delineating blood vessels from surrounding brain tissue enables the quantification of cerebrovascular changes. The PDFs indicate the total blood vessel diameter change in time from day 0 and day 700. Below the circle of Willis, there is a statistical correlation between PDFs and systemic BP (p-value = 0.0007), demonstrating that increased MAP is related to decreased average vessel diameter and PDF. The BP and MAP measurements correlate well with most patients' non-invasive mean PDF diameter measurements below the circle of Willis. Since cerebral changes have been hypothesized to precede systemic hypertension [9,10,26], our methodology may present a tool for potentially initiating early treatment to prevent or optimize management of systemic HBP in conjunction with other approaches including cognitive testing. This finding is important as it suggests that the remodeling of vessels due to increasing blood pressure occurs prior to the onset of diagnosed essential hypertension. Individuals in the current study were explicitly selected to have pre-hypertensive values of blood pressure. Of equal importance, the methodology can determine the relationship of cerebrovascular remodeling to cortical small vessel disease and lacunar lesions known to occur with advanced hypertensive disease and with implications for stroke and dementia [27,28]. The correlation of PDF to MAP was independent of patient gender. The difference in PDF of vessels above circle of Willis was statistically significant with age, which indicates that older patients have constricted cerebral vessels, which may put them at a higher risk of strokes.
In some patients (C and D) the change in CDF did not correlate to changes in MAP, which may indicate impaired auto regulation of cerebral blood flow potentially due to cerebrovascular remodeling [5,26]. The 3D cerebrovascular model reconstruction of these patients demonstrated significant vascular changes between day 0 and day 700. These results may indicate that drug therapies prescribed using systemic BP alone may not provide optimal medical management. Lack of correlation between CDF and MAP may indicate cerebrovascular changes and a higher risk of cerebrovascular adverse events which necessitates more frequent monitoring and/or optimization of medical management, despite having normal systemic BP and MAP. Using a combination of BP and CDF changes may help minimize the occurrence of adverse events.
The proposed approach uses MRA, which directly visualizes overall cerebrovascular anatomy and provides a higher resolution for small blood vessels compared to Doppler Ultrasound which primarily assesses blood flow or the characteristics of a single or small number of vessels within its field of view. MRA is required for visualization of small cerebral blood vessels to obtain an accurate 3D vascular structure. Routine screening using the proposed MRA-based method would be expensive. The clinical application of importance would be treatment of resistant hypertension, which is now poorly understood [29]. A recent study presents convincing data suggesting that such patients may have abnormalities of the vasculature, e.g., of the circle of Willis, that could be detected readily with our technique [30]. Other work has suggested that such hypertension may be due to impingement of the cerebral vasculature on brain stem blood pressure regulatory areas [31]. These abnormalities could be detected with our technique. These vascular bases of hypertension, once identified, would guide clinical treatment in these cases which have not be remedied by currently used pharmacological dosing. Additionally, patients at a higher risk of developing hypertension due to familial or clinical history may be screened using the proposed approach despite the higher cost to detect vascular changes and manage high blood pressure before disease onset.
The segmentation algorithm and metrics for vascular and blood pressure changes (CDF, PDF) are not limited to cerebral vasculature. These methodologies may also be used to quantify vascular changes in other end organs that are sensitive to blood pressure (e.g., kidneys).

Limitations
While our segmentation algorithm significantly improves on automatic segmentation methodologies, it is limited by the resolution limit of the MRI machine performing the MRA scanning. The CDF diameters ( Figure 5) start at 0.5 mm because the distance map calculations determine radius from the edge of a blood vessel and a pixel in the MRA imaging represented 0.25 mm. Any value less than 0.5 mm would not be accurately represented due to the resolution limit. Subsequently, the accuracy of the statistical analysis decreases with decreasing blood vessel size (smaller blood vessels < 10 micrometers) above the circle of Willis).
Various over-the-counter medications and supplements were used by the subjects over the time period of this study; however, the BP changes caused by these medications should be minimal. Nonetheless, larger sample sizes are required to establish definitive relationship with progression to HBP. Despite these limitations, our method is relevant to understanding brain pathology relevant to hypertension whether such pathology precedes or follows the establishment of clinical hypertension.
Although it is hard to ask participants in a study that lasts for more than two years to visit the lab periodically to get their blood pressure measured, we suggest that taking blood pressure readings more frequently (maybe every three months) could provide more accurate observation of the MAP while excluding its possible temporal fluctuations, which could enhance the accuracy of this study.

Conclusions
Alterations in the cerebral vascular tree could be non-invasively detected by the analysis of MRA imaging data. Cerebrovascular changes are correlated to MAP below the circle of Willis. The improved segmentation algorithm coupled with the CDF and PDF estimations could indicate cerebral vascular alterations, which could assist clinicians to propose appropriate medical plans of hypertension.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

I.
Initialization sequentially using EM Algorithm Consider F being the marginal distribution of gray level, EM algorithm [22,32] is used to initially build an LCDG model that approximates F as explained in the following steps: A mixture P K of dominant mode K positive discrete Gaussians (DG) is used to approximate F.
Subordinate components of the LCDG alternatingly approximate the deviations between F and P K as follows: Separating and scaling up the positive and the negative deviations to get the two probability distributions, D p and D n .
Iteratively finding subordinate mixtures of positive and negative DGs using the same EM algorithm. These mixtures should best approximate D p and D n . The mixtures sizes, C p − K and C n are found by sequentially minimizing the error between each distribution (D p or D n ) and its corresponding mixture model with the components number.
Scaling down the positive and negative subordinate mixtures and adding them to the dominant mixture, which gives the initial LCDG model whose size is C = C p + C n .
The initial LCDG has K dominant weights w p,r , where r = 1,2, . . . ,K. The sum of these weights is equal to 1 ∑ K r=1 w p , r = 1 . In addition, there are several lower valued subordinate weights that fulfill ∑ C p r=K+1 w p , r − ∑ C n l=1 w n , l = 0 . II. Refining LCDGs using Modified EM Algorithm The initial LCDG was refined by estimating the local maximum of the log-likelihood in (3) using the DGs adapted EM algorithm in [18], which extends the conventional EM algorithm in [22,32] to handle the alternating components. The described process was shown to be converging to a local log-likelihood maximum, using similar considerations as in [22,32]. Moreover, it was demonstrated in [18] that this process is a block relaxation minimization-maximization.
Considering unit factor constraints in (A2), the log-likelihood in (3) can be equivalently given as: