Normal BOLD Response to a Step CO2 Stimulus After Correction for Partial Volume Averaging

Cerebrovascular reactivity (CVR) is defined as the change in cerebral blood flow induced by a change in a vasoactive stimulus. CVR using BOLD MRI in combination with changes in end-tidal CO2 is a very useful method for assessing vascular performance. In recent years, this technique has benefited from an advanced gas delivery method where end-tidal CO2 can be targeted, measured very precisely, and validated against arterial blood gas sampling (Ito et al., 2008). This has enabled more precise comparison of an individual patient against a normative atlas of healthy subjects. However, expected control ranges for CVR metrics have not been reported in the literature. In this work, we calculate and report the range of control values for the magnitude (mCVR), the steady state amplitude (ssCVR), and the speed (TAU) of the BOLD response to a standard step stimulus, as well as the time delay (TD) as observed in a cohort of 45 healthy controls. These CVR metrics maps were corrected for partial volume averaging for brain tissue types using a linear regression method to enable more accurate quantitation of CVR metrics. In brief, this method uses adjacent voxel CVR metrics in combination with their tissue composition to write the corresponding set of linear equations for estimating CVR metrics of gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). After partial volume correction, mCVR and ssCVR increase as expected in gray matter, respectively, by 25 and 19%, and decrease as expected in white matter by 33 and 13%. In contrast, TAU and TD decrease in gray matter by 33 and 13%. TAU increase in white matter by 24%, but TD surprisingly decreased by 9%. This correction enables more accurate voxel-wise tissue composition providing greater precision when reporting gray and white matter CVR values.


INTRODUCTION
Cerebrovascular reactivity (CVR) is defined as the change in cerebral blood flow produced by a change in a vasoactive stimulus. Methods for assessing changes in blood flow include ultrasound (McDonnell et al., 2013), SPECT, CT (Chen et al., 2006), and MRI (Vesely et al., 2001). Commonly used vasoactive stimuli include inhaled CO 2 (Vesely et al., 2001;Liu et al., 2019), acetazolamide infusion (Settakis et al., 2003), and lowering of blood pressure. One of the more common methods uses inhaled CO 2 during BOLD MR imaging for measurement of CVR. An advance method for delivering quantitative vasodilatory stimuli (Slessarev et al., 2007) has enabled more accurate quantitation of vascular reactivity in healthy subjects and in patients with vascular diseases altering CVR. Despite this, the accuracy of the measurement also depends on post-processing methods. Standard CVR models assume a linear relationship of the BOLD response with the partial pressure of end-tidal CO 2 (PETCO 2 ) for calculating the magnitude of the response. However, the CVR response is more complex than just a scaled version of PETCO 2 . Multiple studies have addressed the problem of modeling CVR responses with development of other methods for taking into account the speed of the response (Poublanc et al., 2015), the time delay (Champagne et al., 2019), and competition from other voxels in the brain .
Another limitation of CVR mapping is the low spatial resolution of the CVR images on 3T MRI. Voxel size can be as small as 2 mm when using a BOLD acquisition with advanced imaging techniques such as multi-band imaging. However, most common voxel size would range from 3 to 4 mm. Such voxels will contain different tissue constituents resulting in partial volume averaging effects. Partial volume effects (PVE) occur when the intensities from gray matter (GM), white matter (WM), and CSF are averaged together within a voxel. However, using the proportion of each tissue constituent within large voxels, it is possible to apply a correction in order to separate the intensities of each tissue constituent. This proportion is derived from tissue segmentation using a high-resolution anatomical image. Partial volume correction (PVC) is critical in patients with cortical atrophy since they are more likely to have fewer pure voxels that contain only one tissue constituent compared to normal individuals. PVC has previously been applied in PET imaging (Rousset et al., 2007), ASL MRI (Asllani et al., 2008;Petr et al., 2011Petr et al., , 2013Binnewijzend et al., 2013), and spectroscopic imaging (Lu et al., 2007;Erlandsson et al., 2012;Bhogal et al., 2020), but not to CVR imaging. Different methods exist to correct for PVE. One of the most straightforward methods is to mask out voxels with tissue fraction below a given threshold. In a study comparing blood flow between young and elderly subjects, voxels with less than 80% gray matter were disregarded (Asllani et al., 2009). However, for a voxel-wise correction, more advanced methods exist such as a Bayesian (Chappell et al., 2011) or a linear regression methods (Asllani et al., 2008). The linear regression method, first introduced by Asllani et al. (2008), is more common and has been used frequently to correct ASL maps.
In this study, we introduce new techniques to further refine the temporal and spatial analyses of CVR images. First, the standard CVR model is compared to the convolution model using the speed of CVR response metric introduced in a previous study. Second, a new model is proposed that utilizes the existing speed metric to measure a time delay (TD). Lastly, a linear regression is applied to correct PVE. This method uses adjacent voxel CVR metrics in combination with their tissue composition to write the corresponding set of linear equations solving for CVR metrics derived from GM, WM, and CSF. This correction is applied to all CVR metrics and results are reported in the major vascular territories of 45 normal subjects before and after correction for PVE.

Subjects
This study was approved by the Research Ethics Board of the University Health Network and conformed to the Declaration of Helsinki. Written informed consent was obtained in all 45 healthy volunteers (20 women). The age range included 20 individuals between 20 and 30 years of age, 20 between 31 and 60 years of age, and 5 over 60 years old. Mean age for all subjects was 39.2 ± 16.5 years. Subjects were in good health, non-smokers, and taking no medication.

Vasodilatory Stimulus
Cerebrovascular reactivity was assessed by measuring the change in BOLD signal to a standardized change in end-tidal partial pressure of carbon dioxide (PETCO 2 ) as the vasodilatory stimulus. PETCO 2 and end-tidal partial pressure of oxygen (PETO 2 ) were targeted independently of each other and of the subjects' minute ventilation and breathing pattern, by using an automated gas blender and sequential gas delivery breathing circuit (Slessarev et al., 2007) (RespirAct;Thornhill Medical, Toronto, ON, Canada). Targeted PETCO 2 and PETO 2 was achieved within 1-2 breaths by administering blends of gasses according to a previously described algorithms (Slessarev et al., 2007). The targeting sequence used in this study was as follows: (1) maintaining resting PETCO 2 for 2 min, (2) a step increase of 10 mmHg PETCO 2 above the resting level for 2 min, and (3) a step down decrease in PETCO 2 for 2 min back to the resting PETCO 2 level.

Image Analysis
MRI and PETCO 2 data were transferred to an independent workstation and preprocessed using AFNI (Cox, 1996;Cox and Hyde, 1997) and Matlab 2015a. Functional BOLD images were volume registered, slice-time corrected, and co-registered to the anatomical T 1 -weighted scan.

Standard CVR Model
PETCO 2 data were re-sampled and time-shifted to the point of coincidence between the rapid changes in PETCO 2 and BOLD signal using Matlab. Magnitude CVR (mCVR) was calculated on a voxel-by-voxel basis from the slope of a linear least-squares fit of the BOLD signal data series to the PETCO 2 values obtained from the RespirAct. For this regression, the rise of PETCO 2 and brain average BOLD signal were matched (Figures 1A1,A2). CVR is expressed as the percent change in BOLD signal per mmHg change in PETCO 2 . BOLD signal drift was controlled by adding a linear trend in time to the model. The mathematical expression of this model is as follow: Variables α × t account for linear BOLD signal drift; b(t) is a constant baseline over time; ε(t) is the residual time series. mCVR is the slope of the regression of the BOLD signal for all the acquired datapoints in each voxel during the stimulus.

Convolution CVR Model
BOLD response was modeled as the PETCO 2 convolved with an exponential decay function exp(−t/TAU), where t is time and TAU is the time constant of the vascular response. TAU was allowed to vary from 2 to 100 s in 2 s increments. TAU derived from the convolved PETCO 2 with the highest Pearson correlation coefficient with the BOLD response was taken as a measure of the speed of response (TAU in Figure 1B1). Steady state CVR (ssCVR) was calculated as the slope of the regression between the BOLD response and the convolved PETCO 2 waveform that provided the highest correlation ( Figure 1B2). TAU is expressed in seconds and ssCVR is expressed as percent MR signal change per mmHg of PETCO 2 .
The mathematical expression of this model is as follows: where α × t, b(t), and ε(t) are the same measures as in (1).
later used in the text and figures (Figures 1B1,B2). Note that ssCVR can be thought of as CVR corrected for the speed of the response TAU. This method is explained in more detail elsewhere (Poublanc et al., 2015).

Time Delay (TD)
The convolved PETCO 2 previously calculated was first shifted 2 s earlier in time to ensure it precedes BOLD signal for all voxels for the brain. It was then shifted to the point of maximum correlation with BOLD response, on a voxel-wise basis (Figures 1C1,C2). This shift represents TD and was allowed to vary between 0 and 5 s, in increments of 0.2 s. The mathematical expression for this is: where α × t, b(t), and ε(t) are the same measures as in (1). Convolved_ PETCO 2 is the modified PETCO 2 found in (2). Due to the small shift in the convolved PETCO 2 , ssCVR value can vary slightly from the value found in (2).

Spatial Normalization and Segmentation
Statistical Parametric Mapping (SPM) 8 software was used to segment and spatially normalize the high resolution T1-weighted anatomical acquisition (Ashburner and Friston, 2005). The segmentation provides probability maps of gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). GM, WM, and CSF masks were generated for each patient by thresholding the corresponding probability maps at 0.7.

Nearest Neighbor Interpolation (NN)
All CVR metrics maps were re-sampled to the T1-weighted resolution using NN interpolation.

Partial Volume Correction (PVC)
The proportion of GM, WM, and CSF within the large BOLD voxels (3.5 × 3.5 × 3.5 mm) in the functional BOLD imaging protocol differs between subjects. Since normal cortical gray matter is approximately 2.5 mm in thickness, all 3.5 mm isotropic voxels over cortical structures (except for a few voxels where gyral surfaces are in contact with each other across a sulcus) will contain less than 100% gray matter (Figure 2). This can introduce an underestimation bias in the voxel content of gray matter, particularly in the setting of cortical atrophy. The following method was applied to correct for this volume averaging. Let's consider the BOLD acquisition voxel V 0 with signal intensity I 0 . Note that I 0 could represent mCVR, TAU, ssCVR, or TD. I 0 is equal to the weighted average of the intensity of GM (I g ), WM (I w ), and CSF (I c ) that corresponds to I 0 = g 0 I g + w 0 I w + c 0 I c, where g 0 , w 0 , and c 0 are the proportion of GM, WM, and CSF inside V 0 . This equation cannot be solved by itself since there are 3 unknowns (I g , I w , and I c ). In order to solve it, surrounding voxels are needed. Under the assumption that the intensity of each tissue constituent is constant within the surrounding voxels, a multiple linear regression model can be used to solve for I g , I w , and I c . In this analysis surrounding voxels were all voxels within a 9 mm radius. This yields 81 voxels including V 0 corresponding to 81 equations: I n = g n I g + w n I w + c n I c + ε n where n = 0,1,... 80.
For each BOLD voxel n, the quantities g n , w n , and c n represent the proportion of GM, WM, and CSF, respectively, in the BOLD voxel. ε n represents the error. The unknown quantities I g , I w , and I c are the intensity of GM, WM, and CSF in the BOLD voxels. The proportions of the 3 tissue classes within the BOLD voxel are obtained from the high resolution T1-weighted segmentation described earlier that has 1 mm isotropic voxels. After solving FIGURE 1 | Diagrams illustrating CVR metrics calculations. (A1) represents one voxel BOLD response (gray) to the step PETCO 2 (red). BOLD signal is regressed against PETCO 2 to obtain its slope, mCVR (A2). The PETCO 2 is then convolved with an exponential function of time characteristic TAU that is chosen to maximize the correlation with the BOLD response (B1). BOLD signal is then regressed against the convolved PETCO 2 to obtained its slope, ssCVR (B2). The convolved PETCO 2 is then shifted to obtain maximum correlation with the convolved PETCO 2 . This shift corresponds to TD (C1) and the adjusted regression is represented in (C2). Notes that for each of the regressions (A2,B2,C2), BOLD signal drift is added to this model as a linear trend. For simplification, it is not represented here.
the linear equations using linear regression, I g , I w , and I c are assigned to the correct tissue class inside V 0 . The same process is then repeated for all voxels of the brain. Note that voxels inside the 9 mm radius but outside the BOLD mask are excluded from the regression. This method was implemented in a script using multiple AFNI functions and is similar to previous work (Asllani et al., 2008;Petr et al., 2013;Ahlgren et al., 2014). Figure 2 illustrates the process. 3dTfitter from AFNI was used for the linear regression. The regression was applied by minimizing the sum of absolute errors, rather than the sum of square since it is less sensitive to outliers. The correction was repeated for each CVR metric. In addition, mCVR was corrected for partial volume averaging using different radii (5, 7, 9, and 11 mm) to illustrate the change in intensity with radius size.

Vascular Territories Template (VTT)
A template of vascular territories, available online and constructed and reported in a 2019 study (Schirmer et al., 2019), was used in this analysis. Those territories included the cerebellum (CER), anterior cerebral arteries (ACA), middle cerebral artery (MCA), and the posterior cerebral artery (PCA) territories. The inverse transformation matrix previously calculated for each subject was used to transform the VTT back to subject original space, rather than transforming the functional maps to MNI. Indeed, warping the functional maps to MNI would introduce re-sampling, leading to a sub-optimal partial volume correction.

CVR Metrics Averages
The average of CVR metrics (mCVR, ssCVR, TAU, and TD) for the PVC and the NN interpolation maps were calculated within GM of each vascular territory, as well as within WM. We attempted to exclude some of the sinovenous structures from the analysis by removing all voxels with mCVR above 0.7%/mmHg. In addition, voxels with negative ssCVR were removed from TAU averages. TAU in those voxels does not represent the speed of vasodilation since those voxels are likely driven by a steal phenomenon from stronger CVR responses in adjacent neighbors. In normal controls, there are only a small number of voxels with this negative response.

Model Selection
Bayesian Information Criterion (BIC) was calculated for the standard CVR model, as well as for the convolution model, For example, voxel V 0 contains 81% gray matter (g 0 ), 19% white matter (w 0 ), and 0% CSF (c 0 ). The intensity I 0 in V 0 is an average of the intensities of GM (I g ), WM (I w ), and CSF (I c ) weighted by their respective matter proportion. The corresponding equation would be I 0 = g 0 I g + w 0 I w + c 0 I c . Under the assumption that each matter intensity (I g , I w , and I c ) is constant within a group of voxels surrounding V 0 , 8 similar equations can be written for V 1 , V 2 , . . ., V 8 to solve for I g , I w , and I c inside V 0 . I g , I w , and I c are then assigned to the corresponding tissue matter inside V 0 (D). Notes that voxels surrounding the central one but outside the BOLD mask are excluded from the regression. In this example, the solution I c is not assigned since there is no CSF inside V 0 . For this example, only 8 voxels surround V 0 . In the actual image analysis, the surrounding voxels were all voxels within a 9 mm radius. in order to compare them. The model with the lowest BIC corresponds to the best model. BIC = n ln(SSE/n) + (p + 1) ln(n) SSE is the residual sum of squares, n is the number of time points, and p is the number of parameters used in the model. The second term in the BIC expression introduces a penalty when adding more parameters to the model. For the standard model (Eq. 1), there are 2 parameters (mCVR and), thus p = 2, and for the convolution model (Eq. 2), there are 3 parameters (ssCVR, TAU, and), thus p = 3. BIC was averaged within gray and white matter separately for each model.

RESULTS
The BOLD response in one voxel of one subject is shown in the diagram of Figure 1. This clearly shows that the rise in BOLD signal following a step hypercapnia stimulus is not instantaneous but rather follows an exponential function. An example of mCVR, ssCVR, TAU, and TD maps generated in one subject are shown in Figure 3 in original space and original resolution (ORIG). The same maps were then resampled onto the anatomical space using nearest neighbor mode (NN) for comparison with the partial volume corrected maps (PVC).
After partial volume correction, mCVR and ssCVR increase in gray matter, respectively, by 25 and 19%, and decrease in white matter by 33 and 13%. In contrast, TAU and TD decrease in gray matter by 33 and 13%. TAU increases in white matter by 24% and TD surprisingly decreases by 9%. Overall, this indicates that PVC effectively draws gray and white matter measures toward their true values. Measurements for gray matter of each separate vascular territory follows the same tendency as for the whole gray matter. Those results are combined in Table 1 and Figure 4. In addition, the results performed on mCVR at different radii sizes were the following: the average GM for NN, PVC with the 5, 7, 9, and 11 mm radii were, respectively, 0.17 ± 0.03, 0.19 ± 0.04, 0.20 ± 0.03, 0.22 ± 0.04 %/mmHg, and 0.22 ± 0.04 %/mmHg. The average WM for NN, PVC with the 5, 7, 9, and 11 mm radii were, respectively, 0.09 ± 0.02, 0.05 ± 0.02, 0.05 ± 0.02, 0.06 ± 0.02 %/mmHg, 0.07 ± 0.02 %/mmHg. Average BIC over gray matter was 76.6 ± 58.7 and 61.2 ± 63.8 for the standard and the convolution model, respectively. For FIGURE 3 | Maps of CVR metrics (mCVR, ssCVR, TAU, and TD) for one subject represented using the original resolution (ORIG) in the first row and processed using the partial volume correction (PVC) in the second row. The raw BOLD image and the high resolution T 1 weighted image are shown in the first column. Using the matter segmentation from the T 1 -weighted high resolution (1 mm isotropic) scan, the original CVR metrics maps of 3.5 mm isotropic resolution (ORIG) were corrected for partial volume effect and mapped onto the 1 mm isotropic grid. 3.0 ± 0.5 2.8 ± 0.9 white matter, average BIC was 58.5.6 ± 50.0 and 48.0 ± 53.1 for the standard and the convolution model, respectively. BIC values for all individual subjects were smaller for the convolution model than for the standard model (p < 0.001).

DISCUSSION
In this study, we provide an algorithm that improves the accuracy of CVR measurements in the brain by reducing volume averaging artifacts arising from mixtures of these three tissue classes within the relatively large voxels used in BOLD MRI. This becomes particularly important in aging individuals and in patients with disease related cerebral atrophy. The degree to which averaging of the different signal intensities of GM, WM, and CFS is controlled reduces misclassification and dependence on voxel size. This correction is particularly important when comparing patients with cerebral atrophy or white matter hyper intensities against a healthy control population (Thomas et al., 2011). The model for this correction assumes uniform intensity of each of the GM, WM, and CSF signal intensities within a small region. However, the smaller the radius, the more unstable the linear regression will be. On the other hand, the assumption of uniform intensity is likely to be false over larger radii. Our results shows that the average GM increases with radius sizes before plateauing at a radius size of 9 mm. The average WM decreases from NN resampling to PVC at 5 mm radius and then increases slightly with radius size. We have chosen a 9 mm radius; however, a simulation would have to be performed to clarify the optimum radius for obtaining the highest precision. In this study, the BOLD signal from CSF was not assumed to be null as it often is when correcting ASL images for partial volume effect. CSF signal was incorporated in the linear regression since large BOLD signal in surface veins can extend to the surrounding CSF. Since this technique requires high resolution T1-weighted images for segmentation, co-registration with the CVR is an important step. Our T1-weighted images were co-registered to the CVR raw images with a rigid body registration and visually inspected. The correction was applied on the CVR metric maps in original space. Vascular territories were transformed from MNI space to original space. This is an important point since applying the correction to the maps in MNI space or transforming the corrected maps in MNI space would have introduced additional spatial re-sampling and smoothing, diverging from the purpose of partial volume correction. The algorithm requires neighbor information that provides a more precise classification accuracy for CVR. Further work is needed to test this algorithm on various disease groups using similar testing parameters, to assess consistency, and to evaluate its impact on improving classification accuracy.
Improving accuracy may be possible by examining the point spread function associated with signal bleeding from the low resolution acquisition. This PSF can then be applied to the anatomical image before segmentation for a better understanding of the low resolution contribution of the different intensities from the three tissue classes. This adjustment can be used before the regression to correct for partial volume effect as shown in spectroscopic imaging in a recent study (Lu et al., 2007). However, considering the size of the BOLD imaging voxels (3.5 mm isotropic) compared to the spectroscopic voxels (5 × 5 × 10 mm), the contribution is expected to have a smaller effect, but it is worthwhile to examine to further improve the current implementation.
Note that the partial volume correction was applied on the parametric images as opposed to the raw BOLD signal. Both approaches should be equivalent when the model is linear such as the one used to generate mCVR. However, we might observe more discrepancy when the model is non-linear such as the one used to extract TAU and ssCVR. Although applying this correction on the raw BOLD signal is processing intensive, this comparison might need to be investigated in the future.
Adjusting for the speed of the response in the model also improves the fit significantly. Indeed, taking the number of parameters used in each model, the BIC for the convolution model is smaller than for the standard model. Although the convolved exponential model is able to fit BOLD responses from healthy as well as diseased vasculature, the redistribution of blood flow in the presence of steal physiology in the setting of steno-occlusive cerebrovascular disease remains problematic. The negative TAU under these conditions no longer accurately indexes the degree of vascular disease. In the setting of steal physiology negative TAU can be paradoxically short. Under these pathological states, the CO 2 stimulus first arrives in the circulation of tissue with intact vasodilatory responses. This circulation lowers its resistance to flow followed by rapid redirection of blood flow away from tissue that cannot change its flow resistance as it is already markedly dilated in attempting to maintain resting blood flow at normal levels. Other methods have used cross correlation time-shift to extract a measure of time delay. However, a cross-correlation technique cannot separate the blood arrival time component from the response time component. In fact, the use of just an amplitude and a shift poorly models the step BOLD response since the main time delay component appears as a slow rise of the response rather than a shift in the response. The solution we provide here is to introduce a time shift parameter (TD) to the convolution model. Note that TD attempts to measure blood arrival time. However, as opposed to ASL and DSC imaging, cerebral blood flow is altered during CO 2 manipulations. Under those circumstances, blood flow will be redistributed, and true blood arrival time might be altered. Therefore, this parameter was not named blood arrival time but time delay. TAU contributes to a much larger component than TD in the model fit. Time delay differences across the brain should not exceed more than 3-4 s (Poublanc et al., 2013), only marginally improving the measurement particularly with a 120 s hypercapnia stimulus. In addition, the solution is limited by the protocol as we only applied a single hypercapnic step, and temporal resolution is low with BOLD images acquired at TR of 2.4 s. The protocol used is better suited to measure TAU rather than TD. Future work may explore the use of multiple short CO 2 steps to improve the performance of TD measurement. In addition, blood arrival time could be altered with increasing hypercapnia since a rise in CO 2 induces an increase in the speed of blood flow. Nevertheless, variations in TD across diseased states relative to healthy control participants allows for group and regional comparisons such as brain matter type, hemispheric asymmetries, and variability across vascular territories. In particular, we propose that this added metric may hold greater promise for assessments pertaining to vascular aging and its correlation to cognitive reserve.
Our rationale for using a standard step PETCO 2 protocol has been guided by the fact that our method for increasing arterial PETCO 2 10 mmHg above resting levels can be achieved in one breath. The sharp step rise and fall in CO 2 enabled clear observation of the BOLD response that led to identification of the exponential rise in signal. Various research groups have also utilized a more attenuated step protocol across numerous published studies (Goode et al., 2009;Anazodo et al., 2016;Leung et al., 2016;Moreton et al., 2016). Over time, the sharp step protocol has, in part, become a standard stimulus largely due to its value in measuring the speed of vasodilation while reducing scan time.
Finally, this report provides reference values of the derived CVR metrics in the vascular territories of the healthy human brain. One limitation is the heterogeneity of age within the study group. However, we have previously shown that, based on our published methods, there are no significant differences in the speed and magnitude of CVR and importantly these values do not differ in GM regions amongst various age groups across the healthy aging continuum, i.e., 18-54 (McKetton et al., 2018). We recognize that there are differences by age documented in some studies (Bhogal et al., 2016;Peng et al., 2018) in the literature. However, the ability to rapidly control and implement precise known changes in arterial CO 2 levels as used in this study attests to the accuracy of the data and the conclusions. However, we note that the sample of participants in our older age group is limited and more investigation is needed to determine the effects of healthy aging over the seventh decade of life.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the University Health Network. The patients/participants provided their written informed consent to participate in this study. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
JP was involved in the development of the scripts development and image analysis. All authors participated in the design and conceptualization of the study, as well as the feedback and writing process following the initial drafting of the manuscript, contributed to the article, and approved the submitted version.