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. 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 final 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% specificity and 99.7% sensitivity in identifying and delineating the brain's vascular tree. The PDFs had a statistically significant 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.

Keywords: magnetic resonance angiography (MRA); blood pressure (BP); systolic pressure; diastolic pressure; arteries; cerebral; hypertension

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], [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], [10]]. 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.

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

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).

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.

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], [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], [21]].

A discrete Gaussian distribution is defined on the set of integers (gray levels) Q = {0,1, ..., Q − 1} by the probability mass function:

$$\psi q|\theta =\begin{array}{l}\mathsf{\Phi }0.5,\hfill \\ \mathsf{\Phi }q+0.5-\mathsf{\Phi }q-0.5\hfill \\ 1-\mathsf{\Phi }\mathrm{q}-0.5,\hfill \end{array},\begin{array}{l}q=0\hfill \\ 1\le q<Q-1\hfill \\ q=Q-1,\hfill \end{array},$$

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

(1) $$p_{w,\Theta }q={{\displaystyle \sum }}_{r=1}^{C_p}{w}_{p,r}\psi q|{\theta }_{p,r}-{{\displaystyle \sum }}_{l=1}^{C_n}{w}_{n,l}\psi (q|{\theta }_{n,l}).$$

The weights w = (wp,1, ..., wp,Cp, wn,1, ..., wn, Cn) are restricted to be all nonnegative and to satisfy

(2) $${{\displaystyle \sum }}_{r=1}^{C_p}{w}_{p,r}-{{\displaystyle \sum }}_{l=1}^{C_n}{w}_{n,l}=1.$$

In general, valid probabilities are nonnegative: $$p_{w,\mathsf{\Theta }}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,\mathsf{\Theta }}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 Fs, its Bayesian estimate F is as follows [[22]]: f(q) = (|R|fs (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:

(3) $$Lw,\mathsf{\Theta }={\displaystyle \sum }\mathrm{q}\in Q\mathrm{f}\mathrm{q}\log p_w{,}_{\mathsf{\Theta }}q.$$

The segmentation algorithm basic steps are following [[16]]:

• (1) For every slice Xs, s=1,.....S,
• (a) First is to gather the marginal empirical probability distribution Fs of gray levels.
• (b) Find a starting LCDG model which is nearing Fs by using the initialization algorithm to approximate the values of Cp−K, Cn, and the parameters w, Θ (weights, means, and variances) of the negative and positive discrete Gaussians (DG).
• (c) Fixing Cp and Cn, refine the LCDG-model with the modified EM algorithm by manipulating the other parameters.(See Appendix A for more details)
• (d) Separate the final LCDG model into K sub models. Each dominant mode has a corresponding sub model. This is done by minimizing the misclassification predicted errors and selecting the LCDG-sub model that has the greatest average value (corresponding to the pixels with highest brightness) to be the model of the wanted vasculature.
• (e) 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.

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 PDF (fX) using $$F_X(x)={\int }_{-\infty }^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.

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.

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).

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).

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], [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]]. 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).

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.

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.

Graph: Figure 1 Framework of the data analysis for quantifying cerebrovascular changes from MRA imaging data.

Graph: Figure 2 (a) Original MRA image slice of sample patient at day 0. (b) Manually segmented ground truth (GT) image from image in (a).

Graph: Figure 3 Illustration of segmentation accuracy and errors of the proposed automatic segmentation (C) by comparing to the ground truth (G) [[5]].

Graph: Figure 4 Example of segmentation algorithm output. (a) Sample image slices of a patient at day 0. (b) Sample 3D reconstruction of the segmented cerebrovascular system using a growing tree model.

Graph: 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.

MAP: Figure 6 Applying a 3D growing algorithm to the volume of binary segmented images allows for visualization of the automatically segmented MRA data. Day 0 of patient C (top left). Day 700 of patient C (top right). Day 0 of patient D (bottom left). Day 700 of patient D (bottom right). The correlation between the CDF and MAP was not found in these patients. Some apparent differences between vascular constructions in areas below the circle of Willis are highlighted.

Table 1 Sensitivity and specificity values for automatically segmented images at day 0, day 700, and cumulative.

Table 2 Statistical evaluation mixed effects linear model, where p-values < 0.05 were considered statistically significant. Diameter denotes size of vasculature in segmentation images. Age, gender, and timepoints are clinically acquired data.

Table 3 BP measurements of patients A, B, C, and D.

Conceptualization, F.T., H.K., A.M., A.S. and A.E.-B.; data curation, F.T.; formal analysis, F.T. and A.E.-B.; investigation, F.T., Y.G. and A.E.-B.; methodology, H.K., A.M., A.S. and A.E.-B.; project administration, A.E.-B.; supervision, S.E.-M. and A.E.-B.; validation, F.T. and S.E.-M.; writing—original draft, F.T., H.K., Y.G., A.M. and A.S.; writing—review & editing, F.T., H.K., Y.G., S.E.-M. and A.E.-B. All authors have read and agreed to the published version of the manuscript.

This research received no external funding.

Data used in this study were collected and approved by the Institutional Review Board of Pittsburgh University.

Informed consent was obtained from all subjects involved in the study.

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

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 PK 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 PK as follows:

Separating and scaling up the positive and the negative deviations to get the two probability distributions, Dp and Dn.

Iteratively finding subordinate mixtures of positive and negative DGs using the same EM algorithm. These mixtures should best approximate Dp and Dn. The mixtures sizes, Cp − K and Cn are found by sequentially minimizing the error between each distribution (Dp or Dn) 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 = Cp + Cn.

The initial LCDG has K dominant weights wp,r, where r = 1,2,...,K. The sum of these weights is equal to 1 $${{\displaystyle \sum }}_{\mathrm{r}=1}^K{w}_p{,}_r=1$$ . In addition, there are several lower valued subordinate weights that fulfill $${\displaystyle \sum }_{r=K+1}^{C_p}{w}_p{,}_r-{\displaystyle \sum }_{l=1}^{C_n{}_{ }}{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.

Consider $$p_{w,\mathsf{\theta }}^{\left[m\right]}\left(q\right)={\displaystyle \sum }_{r-1}^{C_p}{w}_{p,r}^{\left[p\right]}\mathsf{\psi }(\mathrm{q}|{\mathsf{\theta }}_{p,\mathrm{r}}^{\left[\mathrm{m}\right]})-{\displaystyle \sum }_{l=1}^{C_n}{w}_{n,l}^{\left[m\right]}\mathsf{\psi }\left({\mathrm{q}|\mathsf{\theta }}_{n,l}^{\left[m\right]}\right)$$ to be the LCDG at the iteration m. Each signal q ∈ Q contribute relatively to each positive and negative DG at iteration m as specified by the corresponding conditional weights:

$${\mathsf{\pi }}_p^{\left[m\right]}\left(r|q\right)=\frac{w_{p,r}^{\left[m\right]}\mathsf{\psi }(\mathrm{q}|{\mathsf{\theta }}_{\mathrm{p},\mathrm{r}}^{\left[\mathrm{m}\right]})}{p_{w,\Theta }^{\left[m\right]}(q)} ; {\mathsf{\pi }}_{\mathrm{n}}^{\left[\mathrm{m}\right]}\left(l|q\right)=\frac{w_{n,l}^{\left[m\right]}\mathsf{\psi }(\mathrm{q}|{\mathsf{\theta }}_{\mathrm{n},\mathrm{l}}^{\left[\mathrm{m}\right]})}{p_{w,\mathsf{\theta }}^{\left[m\right]}(q)},$$ (A1)

taking into consideration that the following constraints apply:

$${{\displaystyle \sum }}_{r=1}^{C_p}{\mathsf{\pi }}_p^{\left[m\right]}\left(r|q\right)-{{\displaystyle \sum }}_{l=1}^{C_n}{\mathsf{\pi }}_n^{\left[m\right]}\left(l|q\right)=1;q=0,...,Q-1.$$ (A2)

Two main steps are repeated iteratively until the log-likelihood is maximized, the E-step[m] and the M-step[m] which are summarized as follows:

E-step[m]: fixing parameters w[m−1], Θ[m−1], calculate the weights in (A1) from the iteration m − 1.

M-step[m]: maximize L (w, Θ) fixing the weights of (A1) to get the conditional maximum likelihood estimates (MLEs) w[m], Θ[m].

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:

$$\begin{gathered} L(w^{\left[m\right]},{\mathsf{\Theta }}^{\left[m\right]}) = {{\displaystyle \sum }}_{q=0}^{Q}f\left(q\right)\left[{{\displaystyle \sum }}_{r=1}^{C_p}{\mathsf{\pi }}_p^{\left[m\right]}\left(r|q\right)\log p^{\left[m\right]}\left(q\right)- \\ {{\displaystyle \sum }}_{l=1}^{C_n}{\mathsf{\pi }}_n^{\left[m\right]}\left(l|q\right)\log p^{\left[m\right]}\left(q\right)\right]. \end{gathered}$$ (A3)

Using (4), consider replacing log p[m](q) in the first summation with $$\log {\mathrm{w}}_{p,r}^{\left[m\right]}+\log \mathsf{\psi }\left(q{|\mathsf{\theta }}_{ p,r}^{\left[m\right]}\right)-\log {\mathsf{\pi }}_p^{\left[m\right]}\left(r|q\right)$$ and log p[m](q) in the second summation with $$\log w_{n,l}^{\left[m\right]}+\log \mathsf{\psi }\left(q{|\mathsf{\theta }}_{n,l}^{\left[m\right]}\right)-\log {\mathsf{\pi }}_{n }^{\left[m\right]}\left(l|q\right).$$ During the E-step, the conditional Lagrange maximization of the log-likelihood of (A3) under the restrictions of (A2) yields the weights $${\mathsf{\pi }}_p^{\left[m+1\right]}\left(r|q\right)$$ and $${\mathsf{\pi }}_n^{\left[m+1\right]}\left(l|q\right)$$ of (A1) for all r = 1, ..., Cp; l = 1, ... , Cn and q ∈ Q. During the M-step, the conditional Lagrange maximization of the loglikelihood in (A3) under the restriction of (2) and the fixed conditional weights of (A1) results in the DG weights $$w_{p,r}^{\left[m+1\right]}={{\displaystyle \sum }}_{q\in \mathrm{Q}}f\left(q\right){\mathsf{\pi }}_p^{\left[m+1\right]}\left(r|q\right)$$ and $$w_{n,l }^{\left[m+1\right]}={{\displaystyle \sum }}_{q\in \mathrm{Q}}f\left(q\right){\mathsf{\pi }}_n^{\left[m+1\right]}\left(l|q\right)$$ . For each DG, the conventional MLEs of the parameters originating from maximizing the loglikelihood after each difference of the cumulative Gaussians can be approximated with the Gaussian density:

$${\mathsf{\mu }}_{c,r}^{\left[m+1\right]}=\frac{1}{w_{c,r}^{\left[m+1\right]}}{{\displaystyle \sum }}_{q \in \mathrm{Q}}(q · f(q){\mathsf{\mu }}_c^{\left[m+1\right]}(r|q)$$

$${\left({\mathsf{\sigma }}_{c,r}^{\left[m+1\right]}\right)}^2=\frac{1}{w_{c,r}^{\left[m+1\right]}}{{\displaystyle \sum }}_{q \in \mathrm{Q}}{\left(q-{\mathsf{\mu }}_{c,i}^{\left[m+1\right]}\right)}^2·f\left(\mathrm{q}\right){\mathsf{\mu }}_c^{\left[m+1\right]}(r|q)$$

where c can be either p or n. The modified EM-algorithm is true till the weights w become strictly positive. The iterations must be ended if the log-likelihood of (3) becomes almost constant or starts decreasing resulting from rounding errors accumulation.

Associating the subordinate DGs with the dominant terms, the final mixed LCDG-model $$p_C$$ (q) is divided into the K LCDG-submodels $$P_{\left[K\right]}=\left[p\left(q|k\right): q \in Q\right]$$ , one per class k = 1, ..., K, and hence the misclassification rate becomes minimal.