Improved biomechanical metrics of cerebral vasospasm identiﬁed via sensitivity analysis of a 1D cerebral circulation model Biomechanics

Cerebral vasospasm (CVS) is a life-threatening condition that occurs in a large proportion of those affected by subarachnoid haemorrhage and stroke. CVS manifests itself as the progressive narrowing of intracranial arteries. It is usually diagnosed using Doppler ultrasound, which quantiﬁes blood velocity changes in the affected vessels, but has low sensitivity when CVS affects the peripheral vasculature. The aim of this study was to identify alternative biomarkers that could be used to diagnose CVS. We used a 1D modelling approach to describe the properties of pulse waves that propagate through the cardiovascular system, which allowed the effects of different types of vasospasm on waveforms to be characterised at several locations within a simulated cerebral network. A sensitivity analysis empowered by the use of a Gaussian process statistical emulator was used to identify waveform features that may have strong cor- relations with vasospasm. We showed that the minimum rate of velocity change can be much more effective than blood velocity for stratifying typical manifestations of vasospasm and its progression. The results and methodology of this study have the potential not only to improve the diagnosis and monitoring of vasospasm, but also to be used in the diagnosis of many other cardiovascular diseases where car- diovascular waves can be decoded to provide disease characterisation.


Introduction
Among the complications of subarachnoid haemorrhage (SAH) and stroke, cerebral vasospasm (CVS) is the leading cause of cerebral ischemia which leads to death of cerebral tissue cells (Fehnel et al., 2014;Macdonald, 2016). CVS manifests itself in the gradual contraction of intracranial arteries, resulting in a narrowed lumen, and initiates days after SAH. The effects of cerebral ischemia typically become evident after one week. Each year, SAH occurs in between 8 and 10 out of every 100,000 adults, of whom approximately 50% go on to develop CVS, which often results in longlasting neurological impairment or death (British Society of Neurological Surgeons, 2006;Pluta et al., 2009).
The mechanism behind SAH-induced CVS is not completely understood. Several pathological mechanisms have been suggested to explain CVS, including changes in vascular responsiveness, and inflammatory or immunological responses of the vascular wall (Lin et al., 2014). Currently, Transcranial Doppler ultrasound (TCD) is the most commonly used diagnostic tool (Kumar et al., 2016). In patients with CVS, TCD typically shows an increase in blood velocity from baseline values in the affected vessels (Aaslid et al., 1984;Fontana et al., 2015;Harders and Gilsbach, 1987). This biomechanical metric (henceforth referred to as a biomarker) effectively detects CVS in larger arteries, but its sensitivity decreases when CVS affects vessels concealed behind thick bone tissue in the peripheral intracranial vasculature. As a result, CVS is commonly diagnosed by excluding other possible causes and there is a need for more effective diagnostic biomarkers.
Pulse waves in the cardiovascular system originate from periodic contraction of the heart, and propagate through the elastic vessels of the cardiovascular system. These waves reflect at sites of mechanical discontinuity, such as bends, bifurcations, or sudden narrowing, including those caused by vasospasm. Numerical models, in particular 1D and 0D-distributed models, have been https://doi.org/10.1016/j.jbiomech.2019.04.019 0021-9290/Crown Copyright Ó 2019 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). successfully used to capture and describe the physics of propagation of these waves (Alastruey et al., 2007;Reymond et al., 2009).
Several computational models have been proposed to better understand the cause-effect mechanisms governing CVS (Baek et al., 2007;Lodi and Ursino, 1999;Robinson et al., 2010). In these studies, CVS was assumed to affect only the middle cerebral artery (MCA) and the focus was on the currently used biomarker, blood flow velocity, showing its increase in the narrowing vessel. To the authors' knowledge, there is no evidence in the literature of studies that link CVS in different locations (including the peripheral circulation) to pulse waveform features measured proximally.
The construction of patient-specific 1D cardiovascular system models is also challenging because the model inputs (anatomical parameters, boundary conditions, and initial conditions) are difficult or expensive to measure in the clinical setting. It is therefore important to ensure that biomarkers and model predictions are robust under a range of plausible model inputs. An extensive exploration of the model input-space is computationally intensive because the model must be run many times with different sets of inputs, but can be made tractable by replacing the mechanistic numerical model with a fast-running emulator. The most commonly used emulation methods are Gaussian process and polynomial chaos expansions (PCE) (Donders et al., 2015;O'Hagan, 2013). In a previous study, we used a Gaussian process (GP) to emulate a 1D cardiovascular system model to a high level of accuracy . We showed how the GP emulator could reduce the computational cost of model runs by 95%. A GP emulator can treat inputs and outputs explicitly as uncertain quantities, and so by determining the proportion of output variance that could be accounted for by each uncertain input we were able to calculate variance based sensitivity indices for each input and output of the model.
The aims of the present study were to (i) develop a 1D numerical biomechanical model of blood flow in the cerebral circulation, (ii) build a GP emulator of this model, (iii) use the emulator to examine the effect of simulated CVS with varying severity and extension on the properties of the pulse waveform, and (iv) to identify biomarkers that show the greatest sensitivity to the presence and extent of CVS. In the wider context, the present study describes the use of sensitivity indices as a way to identify effective biomarkers, which is a novel approach that has the potential to result in clinically useful tools.

Materials and methods
The study consisted of the following steps: 1. Complete network model: implementation of an arterial network model including vessels typically affected by CVS (complete model). 2. Biomarker-pool selection: identification of pulse waveform descriptors. 3. Sensitivity analysis: development of a Gaussian process emulator using a reduced number of model runs evenly distributed across the input parameter-space. The resulting GP was then used to run a full sensitivity analysis (SA) of the complete model outputs of interest to the input parameters of interest, i.e. vessel radii reduction. 4. CVS simulation: The 1D model was reduced to include only the internal carotid artery (ICA) and its peripheral vessels (reduced model) and employed to simulate CVS with different degrees of severity. The model reduction was done to ensure a constant perfusion of the arteries affected by CVS, and to include the effect of vascular autoregulation. The reduced model included the left ICA, anterior cerebral artery (ACA), and the MCA with its peripheral ramifications. The results of this second set of simulations were then analysed to establish a relation between biomarkers and CVS properties w.r.t. traditional biomarker.

Complete network model
The numerical model used to simulate the propagation and reflection of pulse waveforms in typical vascular networks, together with its validation is described in Melis et al. (2017) and Melis (2017). Briefly, the model is based on the reduced 1D form of the general continuity and Navier-Stokes equations for incompressible flows within narrow straight elastic tubes. The model assumes a flat velocity profile to compute the convective term (which goes to zero) and a parabolic profile for viscous losses. The model outputs are pressure, flow rate, and velocity waveforms. These are time-and space-dependent signals with period T equal to the cardiac cycle (T ¼ 0:8 s). A constitutive equation taking into account the elastic wall behaviour is used to close the equation system whose mechanical properties are set to resemble physiological values for healthy subjects.
The vascular network employed in this study is made of 61 vessels connected in a tree-like configuration (Fig. 1). It starts from the ascending aorta and includes the entire circle of Willis and a detailed description of the MCA daughter branches. The mechanical properties of the arteries were taken from Reymond et al.  Table 1). The circle of Willis network from Alastruey et al. (2007) was extended to include the detailed anatomical description of MCA branches with data from Gibo et al. (1981) and Tanriover et al. (2003).
The numerical solution of the complete network model was computed by means of openBF (Melis, 2018), a finite-volume solver for 1D networks of elastic vessels. Pressure, flow, and velocity waveforms were obtained along five equispaced locations in each vessel. The inlet boundary condition at ascending aorta was set as a typical volumetric flow rate waveform taken from Alastruey et al. (2007) (Fig. 2(a)), and the outlet boundary conditions were set by coupling each 1D model outlet with a three-element Windkessel model of the peripheral vasculature, with typical peripheral resistance and compliance values taken from Alastruey et al. (2007). In order to reduce unrealistic reflections, the first Windkessel resistance value R 1 was changed at each time step to match the 1D segment impedance. The second resistance R 2 was then calculated so that R 1 þR 2 ¼R p , where R p was the peripheral resistance set to a constant value that matched typical resistances Table 1 Cerebral circulation mechanical properties. Posterior communicating artery (PCoA), middle cerebral artery (MCA), anterior cerebral artery (ACA), posterior cerebral artery (PCA), anterior communicating artery (ACoA).

ID
Artery of the distal capillary network at that location. The vessel impedance was then computed as q Â c=A 0 , where q was the constant blood density, c the pulse wave velocity (which is a function of the lumen diameter A ¼ Aðx; tÞ and vessel compliance), and A 0 the constant reference lumen diameter (Table 1).

Biomarker-pool selection
We selected potential biomarkers as waveform features that would be easy to extract in the clinical context, and those that are commonly used to calculate more complex metrics such as the augmentation and pulsatility indices. The biomarker-pool therefore consisted of minimum, maximum, and time average in one cardiac cycle (minðÁÞ; maxðÁÞ, and meanðÁÞ, respectively) of velocity (u) and pressure (P) waveforms, along with their first time derivatives (u t and P t , respectively).

Sensitivity analysis
SA was performed to identify the set of outputs most sensitive to a change of each single input of the numerical model . In a highly non-linear vascular model, waveforms may be sensitive to more than one model parameter or combinations of them. Hence, only outputs highly sensitive to changes in lumen radius were selected as CVS biomarkers. SA was performed considering for each vessel the lumen radius (R 0 ) and its length ('), the wall's Young's modulus (E), and the peripheral resistance and compliance of the Windkessel models at the outlets (R p , and C p ). We computed the sensitivity to R p because this remained constant throughout each model run, whereas the individual Windkessel resistances R 1 and R 2 changed at each time step.
SA was done by means of the analysis of variance decomposition (equation 11 in Saltelli et al. (2010)) from which Sobol's sensitivity indices were computed (Saltelli et al., 2010;Sobol, 2001) using the method described in our earlier work , where the numerical model is replaced by a fast-running GP emulator Oakley and O'Hagan, 2004;O'Hagan, 2006). The GP emulator had a zero mean and a squared exponential kernel as described in our earlier paper , and the GP was fitted using the inputs and outputs from an initial set of 50 numerical simulations as described below.
The sensitivity indices measure the contribution to variance in each model output (i.e., variance of a waveform biomarker) from the variance of each model input (first-order indices), or from the interaction of one input with other inputs (total-order indices). By definition, Sobol's indices are in the ½0; 1 interval; in the case of first-order indices, S i;j ¼ 0 indicates that none of the output j variance can be ascribed to input i, and S i;j ¼ 1 indicates that the output j variance depends only on input i. The interactions between inputs are taken into account by total-order indices, and the difference between total-and first-order indices returns higher-order indices, H j . For example, H i;j > 0 indicates that the output j variance is influenced by the interaction between input i and other input variables. The aim of SA in this case is to identify those model outputs sensitive to changes in radius (high firstorder indices w.r.t. input R 0 ) but not to input interactions. The model outputs satisfying these conditions were selected as CVS biomarkers; the SA was performed in four steps: All model inputs were initialised with typical reference values (Table 1) and then changed within AE50% of their reference values as such variation is plausible for the CVS scenario (Findlay et al., 2016). A homogeneous coverage of the input space was ensured by the use of Latin hypercube sampling method as explained in Melis et al. (2017).

A set of 50 input parameter points was identified by the Latin
hypercube sampling from uniform distributions of each model input in the range AE50% of typical reference values given in (Table 1). A set of numerical simulations run using these inputs were used to train the GP emulator. Latin hypercube sampling ensured an even coverage of the input space, and the number of training runs was based on our earlier experience . 2. The GP emulator was trained using normalised inputs and outputs from the training data, and its hyper-parameters optimised via the stochastic gradient descent method. The number of GP hyper-parameters was twice the number of inputs. 3. The trained GP emulator was used to generate the predictions on the Oðd Â 10 3 Þ (where d is the number of inputs) input set. 4. Sobol's sensitivity indices were computed from these outputs and converted to percent values. The outputs were ranked accordingly to the values of their indices with respect to the radius. In addition, the higher-order indices were computed and used to ensure that only features sensitive to R 0 were selected.
The GP emulator was trained on a set of 50 1D numerical model runs, and was validated against an additional set of 200 model runs. To assess the error made by the GP on predicting numerical outcomes, the mean average prediction error (MAPE) was computed as where N is the total number of input points, and Y i and y i are the i th output from the numerical model and from the GP emulator, respectively.

CVS simulations
To simulate CVS, we used a reduced network model. The reduced model starts from the internal carotid artery and bifurcates into the ACA and the MCA (segments 18, 25, and 23 in Fig. 1, respectively). The inlet boundary condition at the ICA root was set as a typical volumetric flow rate from Alastruey et al. (2007) (Fig. 2(b)). The ACA was included to consider the effect of flow redistribution caused by the CVS and it was closed by a three-element Windkessel model, whereas all the distal vasculature past the MCA was included to simulate the effect of CVS.
The propagation and narrowing levels of vasospasm was decided in consultation with two fully trained Interventional Neuroradiologists from University Hospital of Tours, Tours, France. Various CVS types describing the progression of CVS were considered. These included either CVS spreading from the MCA towards periphery (forward CVS) or from the periphery towards the MCA (backward MCA) in a symmetric or asymmetric configuration (Fig. 3). In all the configurations, the lumen diameter reduction was gradually increased in six stages from the baseline condition (0%) to severe vasospasm (60% of vessel narrowing). Narrowing was uniform across the lumen, with no change in vessel shape.

Complete network model and GP emulator validation
The complete network model was validated against results reported in Alastruey et al. (2007) to ensure correct parameterisation of the additional vessels. The waveforms calculated by the two models match at several locations within the system (Fig. 4a). The complete network model was further validated by comparing predicted velocity values time-averaged over the cardiac cycle in the MCA, and in the presence of vasospasm, with published measurements from SAH patients (Sloan et al., 1989) ( Fig. 4b). The complete network model is in good agreement with literature. Discrepancies are smaller than 7% for lumen diameter reduction < 50%; for narrowing > 50%, the complete network model predictions diverge from the experimental measurements (30% difference at 60% narrowing).
The GP emulator, trained upon 50 simulator runs, was validated against an additional set of 200 simulations. The comparison of numerical results and GP predictions is reported in Fig. 4c where the points scattered along the diagonal line of equality indicate a good agreement between simulations and emulator predictions. The MAPE scored on the test-set was 3:53%, which is small with respect to the 300% whole biomarker variation within the CVS range.
The total simulation time using the complete network model for the 50 point dataset used for emulator training was equal to 31 h on a standard Linux workstation. The emulator training and the SA on the complete d Â 10 3 dataset were performed in 10:4s. The total computational time estimated to perform a complete Monte Carlo analysis using the complete network model would have required approximately 612 h. Therefore, the entire SA by means of GP emulator took 5% of the Monte Carlo computational time.

SA and CVS biomarkers selection
The sensitivity analysis results are reported in terms of Sobol's sensitivity indices (Fig. 5a). The total and first-order sensitivity indices were converted to percentage (i.e. sensitivity indices were multiplied by 100) and reported in bar charts (Fig. 5) as plain and hatched bars, respectively. For instance, the SA indicates that the mean value of the pressure waveform is more sensitive to changes in peripheral resistance (R p ) than in peripheral compliance at the set heart rate (C p in Fig. 5c). The velocity first time-derivative is highly sensitive to changes in lumen radius (R 0 , Fig. 5b); in particular, the minimum value of the velocity time-derivative is exclusively dependent on the changes in radius as the difference between first-and total-order indices is 1%. Therefore, in the fol-lowing analysis, the effectivity of min(u t ) as a CVS biomarker is tested.

CVS simulations
Results are presented as percent changes of the selected CVS biomarkers with respect to the pre-vasospasm value (C%). For different levels of vessel narrowing, the C% was computed as CðxÞ¼100 x where x CVS is the value of the CVS biomarker as the CVS is occurring and x REF is the value of the CVS biomarker in the pre-vasospasm configuration.
The results in terms of time average and minimum first timederivative velocity biomarkers for proximal and distal CVS are reported in Fig. 6. The meanðuÞ biomarker is sensitive to small changes in MCA radius as it increases by 50% for 20% lumen radius decrease (Fig. 6a). The minðu t Þ biomarker is less sensitive to proximal changes in radius and it decreases by 50% only for severe CVS (lumen reduction > 50%). Conversely, in case of peripheral CVS (Fig. 6b), the meanðuÞ biomarker slowly decreases as the lumen radius decrease and a change of 50% is obtained for severe CVS (lumen reduction 60%). The minðu t Þ biomarker change is more sensitive to peripheral CVS as a 50% increase is obtained for 30% lumen reduction in any CVS configuration. In the case of symmetric CVS (circle marker), the increase in minðu t Þ is of 75% for 20% diameter reduction.
In an attempt to identify biomarkers capable of detecting vasospasm from an early stage of inception the velocity biomarkers (meanðuÞ and minðu t Þ; meanðu t Þ, and maxðu t Þ) for 10% lumen reduction for all CVS configurations are reported in Fig. 7. This comparison highlights how different biomarkers can be used for early detection of CVS occurring at different locations. The meanðuÞ and meanðu t Þ biomarkers show a 20% increase in all the cases in which the MCA is affected by CVS (i.e., I; I þ II, and I þ II þ III,i nFig. 7(a, c)), but they only decrease by 8% in case of peripheral CVS (i.e., II þ III and III cases). On the other hand, the minðu t Þ biomarker increases of at least 10% in the peripheral CVS cases (Fig. 7b). The maxðu t Þ biomarker changes are lower than 5% in all the simulated cases (Fig. 7d).

Discussion
The aim of this study was to identify more efficient biomarkers capable of characterising different types of CVS and allow its early detection.
A 1D model of the cerebral circulation was developed to describe the physics of wave propagation through a typical vascular network and to predict the effects of vasospasm on waveforms. A pool of CVS biomarkers was identified through SA on waveform features. The biomarkers examined were pulse waveform features whose variation is caused by changes in mechanical properties due to CVS. A computationally efficient exploration of the parameter space in a 1D model of the cerebral circulation was performed by means of GP emulation.
Our modelling approach showed that a CVS occurring at the measurement location caused an increase (more than 150% in severe cases) of meanðuÞ, as expected, and validated by clinical observations. However, when the CVS occurred more peripherally (Fig. 8), meanðuÞ was only marginally affected by the vessel lumen reduction and was effective only for severe levels of vessel narrow-ing (it decreased by 47% when vessel narrowing was 50%). The decrease was due to the increase in peripheral resistance and, in turn, to the flow diversion towards the distal ICA. This demonstrates that meanðuÞ cannot be a good diagnostic predictor for CVS in distal arteries. Conversely, the minimum gradient of the velocity waveform (minðu t Þ) was more sensitive to CVS affecting distal arteries. An increase (up to 155%) of this biomarker occurred for all the CVS configurations tested. Thus we propose that minðu t Þ has potential diagnostic value for CVS in distal arteries.
The influence of vessel narrowing on the decelerating part of the pulse waveform is what would be expected from simple biomechanical considerations. From the Moens-Korteweg equation (relating wave speed to lumen and other mechanical properties), reduction in lumen diameter will result in a higher wave propagation speed of the reflected waves through the narrowed sections, which in turn will result in an 'earlier' superposition of the reflected wave on the incident waves and therefore on changes to the decelerating slope of the proximally measured waveform. Minimum ðu t Þ occurs immediately after peak systole, when there is superposition of both incident and backward travelling waves.  Figure 1 are compared with waveforms published in Alastruey et al. (2007) and computed in a network without the extended MCA description. (b) Comparison of velocity change in the MCA against percentage lumen diameter reduction for the current study (continuous line) and published experimental measurements (dashed line) (Sloan et al., 1989). (c) Comparison between the GP emulator predictions and the mechanistic model outputs: points lying along the line of equality indicate good agreement between the two methods.  Backward waves are produced by mechanical discontinuities such as a change in radius, and so we would speculate that minðu t Þ responds to changes in the backward waves relative to the incident wave. Further numerical and experimental work would be needed to establish if this is the case.
Overall, our velocity biomarker results were in agreement with previously published experimental work where an increase in meanðuÞ was associated with a decrease in MCA lumen diameter. However, while the numerical model accurately predicted the velocity biomarkers, it lacks a detailed and mechanistic representation of the autoregulation system that is believed to play an important role in the response of the cardiovascular system after haemorrhagic events. In this study, cerebral autoregulation was simulated in an indirect fashion, by imposing a constant flow rate at the inlet of a reduced network. Therefore, a further study should include the entire vascular network (i.e., the whole circle of Willis and both sides of the cerebral circulation) and a dynamic representation of the autoregulation mechanism to ensure a direct estimation of the effects of an extended network on waveforms. This network could also be used to identify better measurement locations for monitoring early CVS onset and progression. Alongside, the numerical results described in the present study can be used as hypothesis to inform experimental measurements, and to validate the newly found biomarker prediction accuracy. Fig. 7. Velocity biomarkers for CVS onset (lumen reduction 10%). Depending on the CVS configuration, different velocity waveform features can be used. CVS including the main MCA branch (first generation, I) causes an increase in meanðuÞ and meanðutÞ, while peripheral CVS (second and third generation, II and III) causes the increase in minðut Þ. Fig. 8. Comparison between the current CVS biomarker (meanðuÞ) and the proposed biomarker (minðut Þ) for the peripheral CVS case. A change in 50% of the minðut Þ biomarker occurs for small changes in lumen diameter (15%) whereas the meanðuÞ biomarker decreases of 50% only in the case of severe vasospasm (60% diameter reduction). The proposed mechanistic model can be used to describe pressure waveforms in transcranial arteries affected by CVS. This model was used to identify more sensitive CVS biomarkers than those currently used in clinical practice. In particular, the mean velocity measured by TCD does not detect distal CVS whereas the minimum gradient of the velocity waveform is capable of differentiating between location and severity of CVS. However, the model used in the present study assumes flat and parabolic velocity profiles for the convective term an viscous losses respectively, and these assumptions may be an important factor that determines whether the CVS biomarkers identified in this study can be used effectively in the clinical setting.
Nevertheless, this study, once clinically validated, has the potential to guide development of better monitoring and diagnostic protocols and technologies for a more effective management of patients not only in the context of vasospasm, but in the treatment of many other cardiovascular diseases, whose presence and progression may be decoded in a similar fashion to provide disease characterisation.

Conflict of interest
The authors declare no conflicts of interest.