Retinal microvasculature and cerebral small vessel disease in the Lothian Birth Cohort 1936 and Mild Stroke Study

Research has suggested that the retinal vasculature may act as a surrogate marker for diseased cerebral vessels. Retinal vascular parameters were measured using Vessel Assessment and Measurement Platform for Images of the Retina (VAMPIRE) software in two cohorts: (i) community-dwelling older subjects of the Lothian Birth Cohort 1936 (n = 603); and (ii) patients with recent minor ischaemic stroke of the Mild Stroke Study (n = 155). Imaging markers of small vessel disease (SVD) (white matter hyperintensities [WMH] on structural MRI, visual scores and volume; perivascular spaces; lacunes and microbleeds), and vascular risk measures were assessed in both cohorts. We assessed associations between retinal and brain measurements using structural equation modelling and regression analysis. In the Lothian Birth Cohort 1936 arteriolar fractal dimension accounted for 4% of the variance in WMH load. In the Mild Stroke Study lower arteriolar fractal dimension was associated with deep WMH scores (odds ratio [OR] 0.53; 95% CI, 0.32–0.87). No other retinal measure was associated with SVD. Reduced fractal dimension, a measure of vascular complexity, is related to SVD imaging features in older people. The results provide some support for the use of the retinal vasculature in the study of brain microvascular disease.

retinovascular signs and individual SVD markers, rather than considering retinovascular abnormalities or SVD features as a whole 12,13,16 .
We aimed to assess an array of retinal measurements in an effort to identify those most closely related to individual and total SVD burden independently of any co-association with common vascular risk factors (VRFs). We hypothesised that associations between retinovascular changes and SVD lesions would be partly explained by, but remain independent of, VRFs. We used quantitative retinal measurements, qualitative (visual ratings) and quantitative indicators of white matter hyperintensities (WMH) (WMH% volume) and SVD burden, medical-history-derived-and measured vascular risk variables, in two cohorts: the Lothian Birth Cohort (LBC1936), and the Mild Stroke Study (MSS). The LBC1936 represents a narrow age range at the 'healthy' end of the SVD spectrum. The MSS represents a wider age range with clinically evident brain vascular disease and a higher burden of SVD and VRFs. Analysis of these separate but complementary population samples allows us to examine the relationship between retinovascular and SVD features at different points on the SVD spectrum.

LBC1936.
In the LBC1936, 866 participants (448, 52% male) of mean age 72.5 years (SD 0.7) attended for clinical assessment; 681 participants completed brain imaging and provided WMH and SVD ratings and volumes. The 603 participants with retinal imaging measurements available for both right and left eyes provided data as the analytic sample for the present study (Table 1). Of the 603 (303, 50% male), 270 (45%) had hypertension, 58 (10%) had diabetes (Table 1), and 31% had moderate or severe periventricular WMH scores with 19% moderate or severe deep WMH scores ( Table 2). Of the 603, 33 participants reported having had a stroke (of whom 11 also had evidence of stroke on imaging), and an additional 51 had imaging-only evidence of a stroke lesion, giving a total with any stroke of 84 (14%).
Prior to structural equation modelling (SEM), initial analyses in the LBC1936 used regression with 16 retinal variables as exposures and 11 brain variables as outcomes. The results of these multiple linear regression, adjusted for age and sex, are presented in Table 3. Three of the 176 associations survived adjustment and false discovery rate (FDR) correction. Decreased arteriolar fractal dimension (decreased branching complexity) was associated with greater WMH volume percentage in both brain tissue volume (BTV) (β = −0.203, p < 0.001), and intracranial volume (ICV) (β = −0.203, p < 0.001), and higher Fazekas scores in periventricular regions (β = −0.157, p < 0.001). The following nominally significant age-and sex-adjusted associations attenuated to non-significance with FDR correction: arteriolar fractal dimension with Fazekas scores in deep white matter regions; venular fractal dimension with WMH % BTV and ICV; venular tortuosity with BTV mm3, Fazekas score in deep white matter, perivascular spaces (PVS) score and SVD score; and arteriolar branching coefficient with deep and superficial atrophy ratings, arteriolar asymmetry factor with deep atrophy ratings, and arteriolar length-diameter ratio with ICV (mm3). The above analyses were re-run after excluding LBC1936 participants with history of stroke (n = 84) (Supplementary Table S1). Two additional associations survived adjustment and FDR correction in this sensitivity analysis; decreased venular fractal dimension was associated with increased WMH % in BTV (β = −0.175, p = 0.03) and ICV (β = −0.178, p = 0.03). The three other previously-significant associations with arteriolar fractal dimension remained significant. Measurement models in SEM to derive the latent traits of VRF, 'WMH load' , 'SVD burden' and retinal 'calibre-complexity' in the LBC1936 sample (Supplementary Figs S1-S3) indicate that almost all measured variables contributed significantly to their hypothesised latent variable (with parameter weights with p < 0.05). The exceptions were: smoking (p = 0.34), as found in the VRF latent variable derived previously in the LBC1936 17 ; and microbleeds in the SVD latent trait (p = 0.29). Nonetheless, these variables were included in the VRF and SVD latent traits. All models had acceptable fit to the data (Supplementary Figs S1-S3).
SEMs used to test hypothesized relationships among latent variables of 'calibre-complexity' and 'WMH load' and 'SVD burden' in LBC1936 are shown in Fig. 1. In addition to the contributions of age, sex and VRFs, 'calibre-complexity' was significantly associated with the latent traits of WMH load (standardised parameter weight −0.181, p < 0.001) (Fig. 1a) and SVD burden (standardised parameter weight −0.141, p = 0.021) (Fig. 1b). This retinal construct of vessel calibre and branching complexity accounted for approximately 3.3% and 2.0% of the variance in WMH load and SVD burden, respectively. However, the association between 'calibre-complexity' and SVD burden attenuated to non-significance once corrected for multiple comparisons (p = 0.054).
SEMs used to test hypothesized relationships among individual retinal variables (identified from age-and sexadjusted regression models, see Table 3), and latent constructs of WMH load and SVD burden in LBC1936 are shown in Fig. 2. In addition to the contributions of age, sex and VRFs, reduced arteriolar fractal dimension was significantly associated with the latent traits of WMH load (standardised parameter weight −0.211, p < 0.001) ( Fig. 2A) and SVD burden (standardised parameter weight −0.168, p = 0.002) (Fig. 2B), accounting for approximately 4.4% and 2.8% of the variance in WMH and SVD respectively. These associations survived FDR correction. The VRF latent trait was not significantly associated with fractal dimension or either SVD or WMH latent trait (Fig. 2).
We tested additional SEMs-to test for co-association of retinal and SVD indicators with VRFs-of the three relationships between retinal variables and individually-measured brain imaging variables that had surviving age-, sex-, and FDR adjustment ( Table 3). The results are shown in Supplementary Fig. S13. All retinal-brain associations remained significant and survived FDR correction.
In age-and sex-adjusted linear regression models, only one association remained significant following correction for FDR; reduced arteriolar fractal dimension was significantly associated with higher WMH scores in deep white matter regions (β = −0.327, p < 0.001). The following 12 nominally significant age-and sex-adjusted www.nature.com/scientificreports www.nature.com/scientificreports/ associations attenuated to non-significance after applying FDR correction: central retinal artery equivalent (CRAE) with WMH % BTV and ICV, variation in arteriolar widths with WMH % BTV and ICV; variation in venular widths with WMH % in BTV and ICV and Fazekas scores in deep white matter; arteriolar fractal dimension with WMH % in BTV, ICV, Fazekas scores in periventricular region, and SVD score; and venular branching coefficient and Fazekas score in deep white matter region (Table 4).
We conducted further multivariable ordinal regression analysis on the surviving association between arteriolar fractal dimension and deep white matter Fazekas score, introducing the VRF component as an additional covariate to age and sex (Table 5). Lower arteriolar fractal dimension remained significantly associated with deep WMH Fazekas scores, with OR (95% CI) of 0.53 (0.32-0.87).

Discussion
To our knowledge, the current study is the first that uses SEM to investigate the relationship between retinal microvascular and brain imaging markers of SVD. Whereas our retinal data reduction approach found most retinal measurements to be largely independent, i.e. not indicators of the same latent trait, we did find four retinal variables with considerable shared variance reflecting aspects of geometry, namely vessel width and complexity. This combined vessel calibre-complexity latent trait was related to SVD burden and WMH load in the LBC1936. Examining retinal variables individually in the LBC1936, we found participants with decreased arteriolar fractal dimension, suggesting of a loss of branching complexity, were more likely to have greater cerebral SVD pathology, independent of VRF. Sparser arteriolar microvasculature was also associated with more WMH. All effect sizes were small. The relatively larger effect sizes of the results using individual retinal variables suggest that the relationship between 'calibre-complexity' and SVD burden is largely driven by the effect of reduced arteriolar fractal dimension. Findings were generally consistent in the MSS participants despite the wider age range and increased burden of vascular disease; MSS participants with sparser arteriolar fractal dimension were more likely to have more WMH in deep white matter. The small non-significant associations between VRF and WMH, accounting for less than 2% of the variance in WMH, are in line with those of Wardlaw et al. (2014) in reporting that VRFs  www.nature.com/scientificreports www.nature.com/scientificreports/ only account for a small proportion of the variance in WMH 17 . Our findings add further weight to the concept that easy to measure, common VRFs explain only a small proportion of both variance in retinal vascular and brain SVD features, leaving a large proportion unexplained.
These findings add further support to the hypothesis that retinal fractal dimension, as a global summary measure of the state or health of the retinal vasculature, might have a significant association with the changes taking place in cerebral small vessels. A loss of branching complexity in regulatory systems of other organs of the body is thought to occur as they become less complex with age 18 . In line with our results, Hilal et al. (2014) found those with reduced arteriolar fractal dimension were more likely to have cerebral microbleeds 12 . Reduced fractal dimension has also been found in Alzheimer's disease [19][20][21] . Decreased branching complexity of the arteriolar and venular networks combined in lacunar stroke, a manifestation of SVD, has previously been reported in the MSS cohort 6 .
A large number of association tests were conducted, and several significant associations became non-significant once corrected for multiple testing. Whereas the amount of associations conducted increases the likelihood of type-1 errors, it is also possible that our robust correction for multiple testing is rejecting modest but true associations; larger sample sizes are needed to confirm them 22 . We interpret all unadjusted significant associations cautiously as they may reflect spurious associations. However, they could provide directions for future research. Without adjustment, some additional retinal measurements, including venular tortuosity and arteriolar measures of branching geometry, appeared relevant for some markers of SVD.
A limitation of the study is its cross-sectional design. Longitudinal studies examining retinal changes and progression of SVD and WMH would be valuable. The smaller sample size of the MSS cohort did not allow for    www.nature.com/scientificreports www.nature.com/scientificreports/ consistent statistical methods across cohorts. Furthermore, despite the use of a validated semi-automated software application; VAMPIRE (Vessel Assessment and Measurement Platform for Images of the Retina, version 3.1, Universities of Edinburgh and Dundee), a degree of subjective human interaction introduces potential for intra-grader variability. Sources of variability of retinal measurements from fundus images have been discussed, among others, in Trucco et al. 23 . We also acknowledge our retinal variables data reduction approach represents a first step towards a composite measure of abnormality and subsequent choice of variables is not definitive. Other studies may reach different conclusions based on their data. The participants' age may have reduced the strength of the associations between VRF, WMH and retinal measurements. The relationship between WMH and hypertension appears stronger in middle age 24 . A systematic review of associations between retinal microvascular changes, dementia and brain imaging features reported stronger findings in middle-aged populations 25 , in line with the age-related attenuation of associations between VRFs, hypertension and cerebrovascular disease and retinopathy 26,27 . Due to the low prevalence in our samples, mircobleeds were included regardless of their location. It is possible that including lobar microbleeds might 'contaminate' this predominantly subcortical SVD score with microbleeds related to cerebral amyloid angiopathy. Future work in other cohorts with greater prevalence allowing for separate analysis by location, lobar versus deep, should seek to determine if microbleeds may play a larger role.
The strengths of the current study include the use of two distinct, yet complementary, well-characterised cohorts. This provides a more comprehensive account of the relationship between the retinal and cerebral microvasculature in different populations, representing different points on the SVD continuum. Furthermore, the LBC1936 as a narrow-age cohort, reduces the troublesome problem where chronological age confounds the influence of other variables. All LBC1936 and MSS participants were imaged using the same retinal camera and MRI scanner, at the same time period, according to the same brain and retinal image analysis methods. This included the use of standardised retinal and cerebral examinations and the use of both quantitative and qualitative WMH measures with the careful exclusion of infarcts from WMH volume. Additional strengths include the use of robust statistics. The sample size in the LBC1936 allowed us to account for measurement error by using SEM to derive error-free latent contracts for VRF, WMH and SVD.
In summary, we provide some evidence to support concomitant retinal and cerebral microvasculature changes, with a small but significant association between reduced retinal arteriolar branching complexity and SVD pathology, known to represent a significant risk factor for dementia and stroke 28 . That sparser vasculature is associated with WMH and SVD burden, suggests that the retinal microvasculature can exhibit corresponding microvascular signs of altered vascular competence leading to tissue damage in the brain. These data extend the research potential of retinal imaging beyond the eye in studying neurovascular disease 1 and paves the way for further studies into mechanism of dementia and stroke and possibly to refine the prognostic potential of detailed retinal microvascular analysis, although we again note that the associations are small. The added value of embedding retinal image analytics in cross-linked health datasets to study SVD is yet to be determined; this requires intra-individual change detection in stratified prospective cohorts to be tracked. As a cross-sectional study, we can only determine associations not causation, nor the progression of cerebral SVD. Despite our efforts in judging the statistical significance of the results with FDR correction, the associations found in this study need future replication in independent samples. Careful validation of our findings is required for better characterizing the retinal changes accompanying SVD, and in considering whether retinal measures have a place in clinical work such as in screening or as outcomes in trials. Further studies are required to elucidate whether changes to retinal properties precede pathology or a response to long-term exposure to risk factors.

Methods
Participants. The LBC1936 29 comprised community-dwelling, mostly healthy older adults, of mean age about 70 years when first recruited in older age. All were born in 1936. The data analysed for the present study, including digital retinal photographs and structural brain imaging, were obtained at a second wave of testing when the participants were approximately 73 years old (N = 866). The recruitment and testing of the LBC1936 has been described in detail previously [29][30][31] .
The Mild Stroke Study (MSS) 5 is a prospective study of patients with recent (within 3 months) clinical lacunar or mild cortical ischaemic stroke. All patients were assessed by an experienced stroke physician. The recruitment, testing and imaging of these patients has been described previously 5,11 .
Both studies were approved by Lothian Research Ethics (LBC1936: REC 07/MRE00/58; MSS: 2002/8/64). The LBC1936 study was also approved by the Scottish Multicentre (MREC/01/0/58) Research Ethics. Written informed consent for participation in both studies was obtained from all participants. The research was carried out in compliance with the Helsinki Declaration.

Measures. Retinal image acquisition and analysis.
In both groups, digital retinal fundus images were captured using the same non-mydriatic camera at 45° field of view (CRDGi; Canon USA, Lake Success, New York, USA). 814 LBC1936 (from Wave 2 of testing) and 190 MSS participants provided retinal images of both eyes. Images were centred approximately on the optic disc. For the present analysis, the retinal images were reanalysed for retinal vascular characteristics using the same semi-automated software package, VAMPIRE by an experienced operator. VAMPIRE image processing and analysis has been described in detail previously [32][33][34] . Briefly, the boundaries of the optic disc (OD) and position of the fovea in a retinal image are first detected and the conventional set of OD-centred circular measurement zones established. Zone B is a ring 0.5 to 1 OD diameters away from the centre, and Zone C is a ring extending from OD boundary to 2 OD diameters away. Next, the software detects the retinal blood vessels present in the image and classifies them as arterioles or venules. The observer, when necessary followed a standardised measurement protocol to preform manual interventions to correct computer-generated labelling of image features, blind to all prior retinal analysis, brain and VRF data. There were www.nature.com/scientificreports www.nature.com/scientificreports/ complete retinal measurements from both eyes for 603 LBC1936 and 155 MSS participants. Rejections were due to poor quality images, eyelashes causing streaks across the photograph, out-of-focus images, and overexposure (in either eye); these occurred in approximately 16% of LBC1936 images, and 8% of MSS images, with an additional 4% of MSS images excluded due to considerable differences in image resolution (arising from deviations from standard operation of the system when imaging).
Sixteen retinal vascular parameters were measured from each image in both cohorts: measures of vessel calibre-central retinal artery equivalent (CRAE), central retinal vein equivalent (CRVE), and the variation in calibre-the standard deviation of arteriolar and venular widths (BSTDa, BSTDv); the gradient of the width of the main arteriolar and venular vessel paths (GRADa, GRADv); measures of branching complexity-arteriolar and venular fractal dimension (FDa, FDv); measures of vessel tortuosity-arteriolar and venular tortuosity (TORTa, TORTv); and measures of arteriolar and venular branching geometry-branching coefficient (BCa, BCv), length-diameter ratio (LDRa, LDRv) and asymmetry factor (AFa, AFv). A lowercase 'a' or 'v' following the variable name indicates a measurement of arteriolar or venular vessels respectively. See Supplementary Material for details on all retinal measurements and how retinal variables were selected for analysis. To reduce the number of variables, reduce multicollinearity and increase reliability, the above-mentioned measurements from both eyes of each participant were averaged to provide a mean measurement for all variables. MRI brain image acquisition and processing. LBC1936 and MSS participants (at time of presentation) underwent brain MRI on the same 1.5-Tesla GE Signa Horizon HDx clinical scanner (General Electric, Milwaukee, WI) with T1-, T2-and T2*-weighted and fluid-attenuated inversion recovery (FLAIR) axial whole-brain imaging. Full details of the brain imaging scanning protocol for the LBC1936 and MSS have been described previously 11,31 . All analyses were performed blinded to all other data. The SVD lesions in both studies were assessed qualitatively and quantitatively using validated methods according to a precursor to the STRIVE criteria 35 . WMH were visually scored using FLAIR-, T1 and T2-weighted sequences on the Fazekas score 36 in both the deep (0-3) and periventricular (0-3) white matter. Appropriate sequences were also rated for the presence of microbleeds (location and number), lacunes (location and number), and perivascular spaces (in basal ganglia and centrum semiovale, 0-4 point score each) according to an established rating protocol 37 . Brain atrophy was coded using a validated template 38 , with superficial and deep atrophy coded separately.
We combined the visual lesion scores into an ordinal 'total SVD score' of 0-4, described previously 39 . Briefly, a scale point was awarded for the presence of (early) confluent deep (2-3) WMH and/or irregular periventricular WMH extended into the deep white matter (3); one or more lacunes; one or more microbleeds; and moderate to severe grading (2-4) of basal ganglia perivascular spaces. These showed face-validity both as an ordinal score and as a latent variable in previous analyses both in the present cohorts and in other studies 39,40 . All rating was performed by a consultant neuroradiologist trained and experienced in SVD features and use of the visual ratings. Quality control of images has been described previously 17,40 .
Quantitative measures of WMH, brain and intracranial volume were obtained using T2*-weighted and FLAIR sequences with a validated semi-automated multispectral image processing tool, MCMxxxVI 41 . This tool was used to measure intracranial volume (ICV, soft tissue structures within the cranial cavity including brain, cerebrospinal fluid, dural and venous sinuses), brain tissue volume (BTV, intracranial volume excluding the ventricular cerebrospinal fluid) and WMH. The structure volumes were measured as absolute values in cubic millimetres (BTV mm3, ICV mm3). Quantitative measures of WMH were expressed as percentage of WMH volume in ICV (WMH % ICV) and percentage of WMH volume in BTV (WMH % BTV).
Covariates. Age and sex were included as covariates in both the LBC1936 and MSS samples. Measures of vascular risk were included as covariates in both samples. VRFs were assessed in the LBC1936 subjects at age ~73 years, at the same session as the retinal photography, and a mean (SD) of 9 (5) weeks prior to brain imaging; they were assessed on presentation in the MSS, at the same time as brain imaging, and approximately four weeks prior to retinal photography. A combination of medical history variables (medically diagnosed hypertension, diabetes, smoking, and hypercholesterolemia), and measured variables (blood pressure [BP], haemoglobin A1c, and plasma cholesterol) were used. The average of three sitting BP measurements were used to derive mean systolic and mean diastolic BP variables in LBC1936 and one BP reading was used for MSS subjects. The above measures were recorded for MSS subjects with the exception of haemoglobin A1c. All measures were performed blinded to all other data. Variables were selected based on a set of measures of vascular risk that we had identified contributed to vascular risk of WMH in previous LBC1936 and MSS analysis 17 .
Statistical analysis. Age-and sex-adjusted linear regression was used to analyse the association between the 16 retinal vascular characteristics and the structural brain imaging-derived measurements in both cohorts. To minimise the potential for type I errors, p values were adjusted according to the false discovery rate (FDR) method 42 . LBC1936 participants with a history of stroke (n = 84, 14%; based on medical history and/or brain imaging appearances) were removed in a sensitivity analysis. Due to the small size and insufficient stroke classification this group could not be divided into stroke subtypes. VRFs were tested as possible explanatory variables for any significant associations between retinal and brain imaging variables, since both retinal vascular abnormalities and SVD features are known to be associated with common VRFs such as hypertension, smoking, diabetes, etc.; this was examined using SEM in LBC1936, and multivariable regression models in the MSS cohort (which we judged to be too small for SEM). See Penke and Deary (2010) for an accessible description of SEM as applied in neuroscience 43 .
The basic questions in the analyses were whether retinal vessel measures were associated with brain imaging measures, and whether these associations were co-associated with VRFs. The LBC1936 is both large in size and has multiple measures of brain white matter health and VRFs. Therefore, in testing the questions above, we were www.nature.com/scientificreports www.nature.com/scientificreports/ able to form multi-variable 'latent traits' (unobservable constructs underlying a combination of correlated individual measured variables) for retinal features, white matter health, and vascular risk. Results from the regression analyses in the LBC1936 were used to motivate the hypothesized relationships subsequently to be tested by SEM.
We showed previously that VRFs, WMH measures and SVD features formed latent variables in the LBC1936 17,40,44 . Therefore, we used the same measurement models to derive the latent variables. Vascular risk was measured as a single latent factor from eight variables; hypertension, diabetes, hypercholesterolemia, smoking, (treated) systolic and diastolic BP, haemoglobin A1c, and plasma cholesterol, as previously 17 . The volume of WMH as a percentage of ICV, and Fazekas ratings in periventricular and deep white matter were used to derive a latent variable of 'WMH load' as previously 44 . 'SVD burden' was measured using a single latent factor with five indicators, namely, Fazekas ratings for both periventricular and deep regions, lacunes, microbleeds, and basal ganglia perivascular spaces, as previously 40 . This was undertaken to test whether including three additional imaging markers of SVD might increase the ability to find significant associations. A single latent 'calibre-complexity' factor was derived from four retinal indicators; two measures of vessel width (CRAE, CRVE), and two measures of branching complexity; arteriolar and venular fractal dimension. The derivation of this latent variable is described fully in the Supplementary Material. All models were estimated using R's lavaan SEM package, version 0.5-22 45 .
Models were estimated using the robust (mean and variance adjusted) weighted least squares (WLSMV) estimator. WLSMV is robust to non-normality and is appropriate for model estimation with categorical data. Standardised regression coefficients (parameter weights, comparable to standardised partial beta weights) were computed for each path in the models. Model fit was assessed using cut-off points of >0.06 for the root mean square error of approximation (RMSEA), and ≥0.90 for the comparative fit index (CFI) and Tucker-Lewis index (TLI). Measurement models for latent traits are shown in Supplementary Figs S1-S3.
We tested the same two questions in the MSS as above for the LBC1936. However, due to the smaller sample size, latent variables were not formed in MSS and we did not use SEM to test hypotheses. Instead, multivariable regression models were applied in the MSS to test for associations between retinal and brain imaging-derived measurements, and the controls were applied for age, sex, and VRFs. To reduce the number of vascular risk parameters and the likelihood of type I errors, principal components analysis (PCA) was applied to the eight measured VRF variables in MSS. The first unrotated principal component accounted for a substantial percentage of the overall variance in VRF variables (26%), with loadings ranging between 0.18 and 0.77, and was used to generate a general VRF score. To validate the use of a principal component score, factor scores for VRFs were derived using the same PCA method in the LBC1936 sample. The correlation between the VRF principal component score from PCA analysis and the VRF latent trait obtained using SEM in the LBC1936 was very strong (r = 0.89). Multivariable ordinal regression analysis was used for WMH and SVD scores in MSS. Results are presented as odds ratio (OR) with 95% confidence interval (CI). Predictors were converted to z-scores, such that the resulting ORs reflect the odds of having higher pathology ratings for each standard unit increase in the predictor variable. Regression analyses were performed with SPSS statistics version 22 (IBM Corp., Armonk, NY).

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.