A precision medicine framework for personalized simulation of hemodynamics in cerebrovascular disease

Abstract

Background

Cerebrovascular disease, in particular stroke, is a major public health challenge. An important biomarker is cerebral hemodynamics. To measure and quantify cerebral hemodynamics, however, only invasive, potentially harmful or time-to-treatment prolonging methods are available.

Results

We present a simulation-based approach which allows calculation of cerebral hemodynamics based on the patient-individual vessel configuration derived from structural vessel imaging. For this, we implemented a framework allowing segmentation and annotation of brain vessels from structural imaging followed by 0-dimensional lumped simulation modeling of cerebral hemodynamics. For annotation, a 3D-graphical user interface was implemented. For 0D-simulation, we used a modified nodal analysis, which was adapted for easy implementation by code. The simulation enables identification of areas vulnerable to stroke and simulation of changes due to different systemic blood pressures. Moreover, sensitivity analysis was implemented allowing the live simulation of changes to simulate procedures and disease progression. Beyond presentation of the framework, we demonstrated in an exploratory analysis in 67 patients that the simulation has a high specificity and low-to-moderate sensitivity to detect perfusion changes in classic perfusion imaging.

Conclusions

The presented precision medicine approach using novel biomarkers has the potential to make the application of harmful and complex perfusion methods obsolete.

Background

Cerebrovascular disease, and in particular stroke, is a major public health challenge. It is a leading cause of death and disability [1]. While there have been advances in prevention and treatment in the past—e.g., mechanical thrombectomy for acute stroke treatment—the overall prevention and treatment results still remain poor [2]. A potential game-changer of stroke treatment success is precision medicine [3, 4]. It aims to provide personalized therapy recommendations based on the individual features of the patient. It utilizes today’s plethora of available patient data as well as mathematical modeling to offer individualized predictions for patients [4]. While highly promising, precision medicine relies on the presence of informative data allowing the differentiation of pathology patterns [4, 5]. In cerebrovascular disease, important information about the severity of stroke risk and potential response to treatment is encoded in individual pathophysiological parameters—in biomarkers—which can be recorded to aid decision-making. Here, one of the most important parameters is the hemodynamic status [6]. This biomarker is already used in a precision medicine approach to identify individual patients benefiting from thrombolysis beyond the currently established treatment time windows which is crucial since often treatment is denied due to time constraints [7]. In chronic cerebrovascular disease, it might aid by identifying areas which are highly vulnerable to stroke [8, 9]. In the clinical setting, however, this data is only available using specialized methodologies, i.e., Dynamic Susceptibility-weighted Contrast-enhanced Magnetic Resonance Imaging (DSC-MRI) perfusion, computed-tomography (CT)-perfusion, arterial spin labeling (ASL) perfusion or functional MRI [10,11,12,13,14]. These techniques may harm patients through contrast agents, significantly prolong the time to treatment and lead to increased costs. Also, standardization of these complex methods is highly challenging [6, 15, 16].

An alternative approach to derive biomarkers for precision medicine is the transformation of routinely acquired data by mechanistic simulations [17]. These simulations integrate domain knowledge by mathematically describing known disease-driving core processes [17]. Interestingly for cerebrovascular disease, several works in the past have developed general mechanistic simulations of the blood flow in the brain [18, 19]. These simulations have the potential to become a contrast agent-free biomarker of hemodynamics for the diagnosis and treatment of cerebrovascular diseases. However, for these simulations, personalization on an individual patient level is still pending making it not applicable in a clinical setting.

Thus, the novel idea presented in this work is a software framework to transform routine structural vessel imaging data as an input to a mechanistic simulation of individual hemodynamics for a given patient. The unique vessel configuration of each patient can be used to simulate hemodynamics to potentially identify areas that are vulnerable in case of stenosis and occlusion. Several use cases can be envisioned for such a framework. It could allow assessment of stroke risk, pre-operative simulation of interventional success like thrombectomy in acute stroke, preventive or therapeutic endarterectomy and stenting of brain-supplying vessels, respectively. Another highly interesting, if rather rare case is the simulation of the outcome of extracranial–intracranial (EC–IC) bypass surgery, e.g., in Moya-Moya disease. Here, there is a special need to predict the success of the surgery [20]. Lastly, the simulation information could be used for the prediction of stroke outcome in conjunction with other clinical and imaging parameters enabling clinicians with an objective criterion for decision support in the acute setting.

Thus, the objective of the presented work was to provide a framework allowing the incorporation of individual structural vessel data to simulate areas of higher hemodynamic vulnerability as a disease biomarker. For this purpose, we developed a pipeline consisting of the following sequential steps: (1) segmentation of vessel information from structural data, in our case from time-of-flight (TOF) magnetic resonance imaging (MRI). (2) Annotation of the vessel tree with an easy-to-use graphical user interface (GUI). And (3) simulation where results can be inspected, and different blood pressure scenarios can be simulated by the user.

The simulation was implemented as a steady-state zero-dimensional lumped model of the Circle of Willis (CoW) and major brain artery circulation. We included individual vessel resistances by 1-dimensional calculation using the individual length and the width of the arteries from patient structural vessel imaging. Hemodynamic measures were calculated using an adapted version of the modified nodal analysis (MNA) [21] coined AMNA, where we simplify the solution of the matrix equations facilitating easier implementation and faster runtime. We show the implemented framework in detail—including a video of the annotation process—and perform an exploratory visual analysis to compare simulation results with perfusion imaging in 67 patients with cerebrovascular disease.

Results

The framework architecture

We developed a novel framework for enabling processing of routine imaging (DICOM) into a simulation tool that provides additional information about the state of the hemodynamics of an individual patient.

Figure 1 gives an overview of the subsequent steps of the framework:

• 1. A structural DICOM TOF-image (Fig. 1a) is processed via segmentation into a binary vessel mask (Fig. 1b).

• 2. This image is then skeletonized (Fig. 1c).

• 3. The vessel skeleton is imported in our annotation tool, in which the vessel tree is displayed in a 3-dimensional fashion and

• 4. Vessel annotation is performed (Fig. 1d and e).

• 5. The annotated vasculature is then transferred to the simulation tool (Fig. 1f).

Fig. 1

The input is a binary vessel (b) mask segmented from a structural TOF-image (a). This image is then skeletonized (c). This skeleton can be imported in our annotation tool, where the vessel tree is displayed in a 3D fashion and annotation is performed (d and e; red = left, blue = right, purple = median). The annotated vasculature is transferred to the simulation tool (f)

The annotation module

The annotation tool is made up of two main components. The segment annotation area, where the vessel segments can be chosen and the 3D view, where the imported vessel tree can be manipulated. We implemented 22 segments of brain-supplying arteries, namely 3 for the carotid artery (common, internal and CoW segment), 7 for the MCA (M1, M2 superior, M2 inferior, M3 superior superior, M3 superior inferior, M3 inferior superior, M3 inferior inferior), 5 for the ACA (A1, A2, A3 inferior and A3 superior), 4 for the PCA (P1, P2 and P3 inferior and P3 superior), as well as the basilar artery, the vertebral artery and the anterior and posterior communicating arteries. There is the option to add bypass vessels or collateral vessels manually, e.g., for planning of interventions or surgical procedures.

Next to segment also the following additional labels can be chosen: “pre-occlusion”, “post-occlusion”, and “occlusion” (in case pre- or post-occlusion cannot be determined with certainty). To record occlusion next to the type of vessel is important for the following simulation step, as subsequently the simulation will ignore segments with this flag. A video footage of the annotation process was uploaded to zenodo [22].

The simulation module

The simulation itself consists of a graphical user interface, that is divided functionally in the toolbar (Fig. 2b I.), the simulation area (Fig. 2b II.), the pressure selection area (4B III.) and the view selection area (Fig. 2b IV.). In the toolbar, simulations can be loaded as well as the type of resistance calculation (here, the default is the resistance calculation presented in the methods section). In 4B III., the pressure boundary conditions can be chosen, for one the blood pressure which determines the driving pressure of the whole system.

Fig. 2

The simulation GUI. It is divided into the toolbar (I.), the simulation area (II.), the pressure selection area (III.) and the view selection area (IV.). As shown on the right-hand side (c), additional information can be extracted for each edge (a)—in this case a left ICA—and each supply area (b)

In addition, also the intracranial pressure, which is kept constant for our use case at physiological parameters, but can be increased to simulate conditions with increased intracranial pressure, e.g., in hemorrhagic stroke. In the view selector area, in 4B IV., the normal view and the pressure view can be chosen (please see examples in Fig. 3). Finally, in the simulation area, the individualized simulation of the vasculature can be inspected. Via right click it is possible to derive more information about an edge or a node (Fig. 2a and c, respectively). Supply areas are color coded according to the calculated pressure (all values in mmHg: > 70 pink, < 70 yellow, < 60 orange, < 50 red). Below 50 mmHg we consider an area being vulnerable to ischemia due to the limits of cerebral auto-regulation (see Eq. 5).

Fig. 3

In our visual comparison we found promising initial results for the detection of vulnerabilities by mechanistic simulations (see Table 1). In this figure, we exemplify this in two patients. Simulation results are shown on the left and DSC-imaging results are shown on the right. In the patient with a proximal common carotid occlusion on the right (upper box), no hemodynamic changes are visible in the DSC perfusion imaging maps indicative of a sufficient CoW to provide all brain areas with sufficient perfusion. The simulation shows correspondingly no areas reaching critical limits, at both normal and lowered blood pressure levels. In the patient with the right MCA stenosis (lower box), the DSC perfusion maps show delay in TPP, MTT and Tmax maps as well as a lower CBF in the right MCA area. This is in line with our simulation which shows critical pressure levels in the M2 inferior node and lowered pressure already at a normal blood pressure of 120/80. When blood pressure is lowered to 90/60, both M2 supply areas become critical in contrast to the other branches. Supply areas are color coded according to the calculated pressure (all values in mmHg: > 70 pink, < 70 yellow, < 60 orange, < 50 red)

Table 1 Detailed results of the visual vulnerability analysis

Comparison of DSC perfusion imaging and simulation results

67 patients (mean age 57 years) were included. 68% had previous stroke and 18% had a previous transitory ischemic attack (TIA). The visual analysis revealed a very high specificity and a moderate-to-low sensitivity for the simulation to detect vulnerability defined by visual DSC-rating. Results averaged over both hemispheres for each perfusion territory were (sensitivity/specificity): ACA 0.42/0.94, MCA 0.47/0.98, PCA 0.16/0.94. For detailed results please see Table 1.

Sensitivity analysis

An example of the sensitivity analysis is shown in Fig. 4, where we give an example how the still normal hemodynamics in a patient with a missing ICA would be affected by additional vessel stenosis. The sensitivity analysis successfully allows individualized simulations of potential disease progressions and blood pressure scenarios.

Fig. 4

Illustration of the sensitivity analysis. The patient has a right-sided ICA occlusion. All simulation results are calculated under normo-tension (120/80 mmHg). In the current state (a), the patient has no changes in hemodynamics. Adding a 90% stenosis of the ICA on the other side (b, red arrow) leads to only a little drop in perfusion pressure in both MCA supply areas (yellow supply areas). This can be attributed to a very marked posterior communicating artery on the left side. Adding a 50% stenosis of the left Pcom leads to high vulnerability of the anterior circulation system (red supply areas) of both sides (c, orange arrow). The color changes within each segment are graphical features indicating the natural flow direction of the blood and are not related to simulation parameters

Discussion

We present the first comprehensive precision medicine pipeline for cerebrovascular disease that is capable to process routine stroke imaging DICOM images for simulation of various boundary conditions to identify brain areas vulnerable for ischemia.

Our framework provides annotation of the arterial vasculature derived from neuroimaging followed by zero-dimensional individualized simulation of brain hemodynamics. Implementation was performed by decreasing the computational burden by a modified MNA, the development of an easy-to-use web user interface-frontend for annotation and a java-based cross-platform simulation tool. An exploratory validation analysis comparing our simulation results with DSC perfusion in patients with steno-occlusive disease revealed promising initial results for the simulation-based detection of vulnerable areas. Our results suggest that mechanistic simulation of blood flow derived from routine structural imaging can serve as an individual biomarker for patients with cerebrovascular disease and might be an alternative to complex and potentially harmful perfusion techniques.

Stroke is a complex disease with a dynamic progression. The initial infarct area—characterized by rapid neuronal loss—is called the core, which is surrounded by tissue that is slowly surrendering to ischemia but is still salvageable. The latter area is coined the penumbra and defines the therapeutic target in acute stroke management [6]. For the understanding of stroke and its implications on treatment strategy, it is essential that the speed of the penumbra-to-core transformation varies greatly and is highly individual. In particular, for weighing benefit and risk for stroke treatment the high interindividual variance of brain cell death is crucial: some patients do not have salvageable brain tissue already a few hours after stroke, whereas in others penumbral tissue was found up to 17 h after stroke [23]. This highly individual and variant stroke progression stands in stark contrast to the current “one-size-fits-all” treatment approach in stroke where patients receive treatment based on guidelines usually only within a predefined time window of up to 4.5 h after stroke [6, 24]. These time windows were established by a statistical benefit-to-risk calculation after lumping together all stroke patients with assumed common pathophysiology.

While it is true that there is a net profit for patients when treatment is applied within this time window, it is also obvious that many patients do not receive treatment who would benefit from it and at the same time patients receive treatment subjecting them to risk of intervention such as bleeding without its benefit. This is due to the above-mentioned fact that, in reality, several stroke subpopulations exist. Here, precision medicine accounts for the individual features and will improve outcome by personalizing treatment [25]. Precision medicine utilizes mathematical techniques and available digital data to provide individualized predictions for patients [3]. It relies on the presence of informative data allowing the differentiation of pathology patterns [4]. In stroke, it has been shown that measuring the penumbra through perfusion as a surrogate is one of the most promising approaches [6, 7]. And indeed, perfusion imaging-based selection for treatment beyond the established time windows is an evidenced precision medicine approach in stroke [7]. For the selection process, predictive modelling might be also applicable [5]. A drawback of this approach, however, is the application of perfusion measurement techniques which are potentially harmful through contrast agents, inevitably prolong the imaging time and are problematic to standardize across centers [6, 15]. Similar considerations apply to chronic steno-occlusive disease. These are patients with continuously worsening symptoms of atherosclerosis who have a high likelihood for a future stroke event. In these patients, potentially harmful perfusion imaging techniques should not be used a priori. Contrast agent-free perfusion imaging methods can be used but are—as mentioned above—hard to standardize. Thus, alternatives are warranted, and mechanistic simulations are promising methods. Here, the relevant (patho)physiological biomarker is not directly measured, but mathematically inferred from conditions recorded through other measurements.

As suggested by our work, in the case of cerebrovascular disease, we can infer information about hemodynamics from the individual vasculature of a given patient. We successfully built a pipeline that can extract the vessel information by segmentation, allows annotation of the vessels and simulates hemodynamic information which we were able to relate to clinical DSC-perfusion imaging through an exploratory visual comparison of DSC perfusion and simulation results. While these results need gold standard validation, they pave the way for further development of techniques that might make the need for perfusion imaging in cerebrovascular disease obsolete for some patients while still providing the necessary information for precision medicine selection of patients for prevention and treatment. Our results are in line with other recent exploratory validation studies [26].

Given that mechanistic simulations work on a priori assumptions about the biological system and perfusion measurements actually record dynamic information, it is unlikely that the information provided by both systems will always be a complete match. This is also evidenced by the low sensitivity and the high specificity. DSC-MRI is sensitive to very small changes in perfusion, whereas a mechanistic simulation is expected to distinguish between relevant categories. Thus, the simulation was not able to pick up on every change noticed by the readers (low sensitivity), but where the simulation found vulnerabilities they were almost always accompanied by corresponding changes in DSC-MRI (high specificity). While the clinical relevance needs to be validated in further studies, our results suggest, that the information might be intersecting enough to allow treatment-relevant predictions on an individualized patient level, and consequently to avoid harmful imaging procedures. Thus, when the question is for example about a general status, i.e., “is there a general vulnerability in the right MCA area for ischemia”, mechanistic simulations might be able to provide this information instead of direct perfusion measurements. Also, since many variants of the CoW exist, the simulation could allow the identification of patients with high-risk for stroke owing to their individual CoW configuration [27].

Another potentially big advantage of mechanistic simulations is the possibility to simulate interventions. In chronic steno-occlusive disease, like carotid stenosis or Moya-Moya disease, potential lumen reopening interventions or EC–IC bypass surgery can be performed. With our solution as presented in this work, it would be feasible to simulate the reopening of a vessel and thus simulate the post-intervention status.

With our framework, it is possible to simulate the response of the vasculature to changes in blood pressure. This can potentially be highly important not only for the determination of areas-at-risk for ischemia, but also to predict the response to interventions and surgery, e.g., blood pressure drops during surgery. This is not possible with perfusion measurements, which can only provide a snap-shot of the status quo. Approximations can be done with acetazolamide challenge measurement [28, 29], but this requires repeated measurements, the application of a drug, and can only be performed within the physiologically tolerable range. A clear advantage of direct dynamic perfusion measurements, on the other hand, is with high likelihood still the recording of subtle changes and small lesions. Importantly, we thus do not claim that mechanistic modeling might be able to make all perfusion measurements obsolete.

As a limitation of our work, our results are exploratory and hypothesis-generating [30]. While our exploratory validation yielded promising results, further validations are needed. We believe, however, that our results are motivating to boost the translation of the work done in the past on the translational development of mechanistic modeling of hemodynamics into the clinical setting. We implemented a 0-dimensional (D) model of hemodynamics which exploits the similarities of such a network to an electric circuit. Next to these 1D and 3D models exist; for an overview of existing methods see Leguy et al. and Perera et al. [31, 32]. While 1D and 3D models are more suited to model local changes, 0D models are more suited to model the general vasculature, but they can be combined to provide complementary information [33]. Here, our framework builds on existing work, but adds (A) the calculation of the resistance over segments with variable diameters. (B) an easy to code-adjusted implementation of the modified nodal analysis which reduces computational demands and (C) a graphical user interface tailored for inspection and manipulation of the simulation. This facilitates the application of such a framework, which is promising as there is much promise in mechanistic modeling of blood flow and perfusion for clinical applications in cerebrovascular disease. There is, however, a need to personalize these approaches.

Our work has other limitations. First, the simulation values were not compared to a gold standard. However, it is very difficult to derive individualized gold standard values for arterial flow. Other noninvasive methods are either non-gold standards themselves, like MRI flow measurement methods, and/or cannot access the complete vasculature, like Doppler-sonography. Intra-operative direct vessel flow measurements recorded during EC–IC bypass surgery might be an option for gold standard measurements—which we will explore in the future—but were not available for the current study. Second, we would like to point out that all relevant pre-processing steps in the pipeline—segmentation, skeletonization, and annotation—were done manually in the pipeline. With the advent of powerful machine learning segmentation methods in recent years, it is very likely that these steps can be automated with sufficient performance. For segmentation, our group has just recently presented deep learning methods to segment the vasculature from structural scans with very high accuracy [25, 34]. The application of deep learning for skeletonization and automated annotation is a current focus of our group. We are thus confident that for future potential applications in the clinical setting simulation results will be obtainable in real-time, at the scanner console, in a few years.

Conclusion

We present the first precision medicine pipeline for cerebrovascular disease that allows annotation of the arterial vasculature derived from structural vessel imaging followed by personalized simulation of brain hemodynamics. This enables further development of precision medicine in stroke using novel biomarkers and might make the application of harmful and complex perfusion methods obsolete for certain use cases.

Methods

Data accessability

The datasets presented in this article are not readily available because data protection laws prohibit sharing the imaging data used in this study at the current time point. Requests to access the datasets should be directed to ethikkommission@charite.de.

Patients

Sixty-seven patients with steno-occlusive disease from an imaging study of cerebral perfusion in stroke patients (PEGASUS study [8, 9]) were evaluated as examples for the presented framework. Patient characteristics can be found in Table 2. This study was approved by the institutional ethics committee of Charité Universitätsmedizin Berlin and the patients gave written informed consent.

Table 2 Patient characteristics

The framework pipeline

The framework pipeline consists of these chronological steps: (a) segmentation, (b) skeletonization, (c) annotation, and (d) simulation. In the following, each of the steps is described in detail.

Segmentation

Structural MRI imaging consisted of time-of-flight (TOF)-MR images. To save time, all images were first pre-segmented using a high-performing published neural net segmentation model [25]. The pre-segmentation model was applied without modifications. All information about the data employed in the training, the architecture, the training regime and the performance can be found in the open-access publication of Livne et al. [25]. The segmentations were then manually corrected to derive ground-truths standards. The segmentations are voxel-based binary representations of the vessel tree for each individual patient. The imaging parameters for TOF-MR-images in the PEGASUS study were: voxel size = (0.5 × 0.5 × 0.7) mm [3]; matrix size: 312 × 384 × 127; TR/TE = 22 ms/3.86 ms; time of acquisition: 3:50 min, flip angle = 18 degrees.

Skeletonization

The segmentations were skeletonized using the DtfSkeletonization module (DtfSkeletonization—MeVisLab documentation) of MeVisLab (website: MeVisLab) [35]. Here, a one voxel skeleton of the vessel midpoints is created with the radius encoded in the voxel value. This skeleton volume is then transferred to the manual annotation module.

Annotation

The annotation module was developed and coded from scratch by the in-house development team.

Within the annotation framework, the transferred skeleton volume is transformed into a Java-based tree structure—the so-called skeleton graph—representing the skeleton as a set of edges containing all necessary geometric information for 3D rendering, and a set of junction vertices. This tree structure is loaded via a RESTful service interface in JSON format into a JavaScript-based web frontend where it is rendered by using the Three.js library as a rotatable and zoomable 3D view.

Within the 3D view, it is possible to select edges and tag them with an item from a list of vessels (artery) descriptors and their anatomic location (visualized by color) interactively. The triple consisting of an edge, a tagged vessel item, and an anatomic location is defined as an annotation. A set of made annotations can be saved in the backend via a RESTful service interface within an appropriate Java presentation. In the final step, based on a skeleton graph and a belonging annotation set, a simulation model serving as the input for our simulation component can be created once the annotation is finished.

Within our framework, it is possible to annotate 3rd order branches (A3, M3, P3). For exemplary patients in our study, annotations were performed until 2nd order branches (A2, M2 and P2), since higher order vessels are unlikely to play a crucial role in steno-occlusive disease.

Simulation

The simulation was developed and coded from scratch in-house. Our model describes the cerebral vascular tree by a planar graph, which is given in Fig. 5. The blood flow through the vessel tree is modeled in analogy to electric circuits in a modified nodal analysis.

Fig. 5

The circle of Willis represented by edges, nodes and supply areas. The alphanumeric abbreviations and AcomA and PcomA stand for certain brain vessels (see legend) represented by edges. The round shapes represent nodes, whereas the quadratic shapes represent supply areas. M1: A1: anterior cerebral artery (ACA) segment 1, A2: ACA segment 2, M1: middle cerebral artery (MCA) segment 1, P1: posterior cerebral artery(PCA) segment 1, P2: PCA segment 2, PcomA and AcomA: posterior and anterior communicating arteries

Edges in the graph represent blood vessels, while nodes represent either supply areas, blood sources or junctions between nodes. Supply areas are marked with square-shaped nodes, initial (round) nodes are blood sources and the rest of the (round) nodes are junctions. The arrows of the edges indicate the flow directions to the supply areas or away from source nodes. The blood flow into the supply areas is provided by the outgoing segments (A2, P2 and M1 or M2) of the circle of Willis.

Modeling the flow and vessel network

The cerebral vessel tree retains an overall Reynolds number allowing to describe the cerebral blood flow in terms of a Newtonian fluid. The arteries are modeled as perfect cylinders and the decrease in blood pressure $\Delta P$ along a cerebral artery of length L, with radius r, blood dynamic viscosity μ for a volumetric flow rate Q is determined according to Hagen–Poiseuille equation:

$$ \Delta P = \frac{8 \cdot \mu \cdot L \cdot Q}{{\pi r^{2} }}. $$

(1)

Hagen–Poiseuille equation is equivalent to Ohm’s law. Therefore, the resistance of an arterial vessel can be defined as:

$$ R = \frac{8 \cdot \mu \cdot L}{{\pi r^{2} }}. $$

(2)

In the presented use case, the vessel diameter is not homogeneous across the whole vessel segment and consequently Eq. 2 is invalid. To account for variable diameters over a segment, the fluid dynamics is applied on infinitely small segments with a constant radius to yield the equation for non-constant radii. The resistance can be therefore derived using the following integral equation:

$$ R = \int_{L} {\frac{8\mu }{{\pi r(l)^{4} }}} dl. $$

(3)

In practice, discrete application is used, as a segment is defined by a diameters-vector of length n−1 = the number of voxels in the segment. To account for the fact that an antiderivative can only be determined for segments that can be constantly differentiated, we need to approximate the condition of continuity. For this purpose, the radius $r(l)$ is described by a continuous linear extrapolation function, which connects the radius of a given voxel with the next voxel in a linear fashion to allow the calculation of (3):

$$ R_{ext} = \sum\limits_{i = 1}^{n - 1} {\left\{ {\begin{array}{*{20}l} { - \frac{8\mu }{{3\pi }}(\frac{{r_{i + 1} - r_{i} }}{{l_{i + 1} - l_{i} }})^{ - 1} (r_{i + 1}^{ - 3} - r_{i}^{ - 3} ){\mkern 1mu} ,\;r_{i + 1} \ne r_{i} } \hfill \\ {\frac{8\mu }{\pi }(l_{i + 1} - l_{i} )(r_{i} )^{ - 4} {\mkern 1mu} ,\;r_{i + 1} = r_{i} } \hfill \\ \end{array} } \right.} {\mkern 1mu} ,\;n \in {\mathbb{N}}, $$

(4)

where $R_{ext}$ stands for the extrapolated resistance.

According to the mass flow law, the amount of blood entering a node must equal the amount of blood that leaves a node (see illustration in Fig. 6a). Mathematically this is described according to Kirchhoff’s first law as follows:

$$ \sum\limits_{i = 1}^{N} {Q_{inp,i} } = \sum\limits_{j = 1}^{M} {Q_{out,j} } , $$

(5)

where $Q_{inp,i}$ is the input blood flow from source i to node n, for N input sources, and $Q_{out,j}$ is the output from node n through vessel j, for M outputs.

Fig. 6

a Graphical illustration of Kirchhoff’s first law. The current entering a node = i₂ + i₃ must equal the current leaving that node = i1 + i4. See blood flow analogy Eq. 5 for the mass flow law. In the figure i correspond to current, R represents resistance and v stands for the voltage. b Illustration of a normal node representing a vessel junction. Each numbered circle in the figure represents a node and each arrow represents a directed edge. The figure illustrates a normal node—numbered as 3—with two incident edges (G13 and G23) and three outgoing edges (G34, G35 and G36). The resulted G matrix row vector for node 3 is: $G(3,3) = G_{3,1} + G_{3,2} + G_{3,4} + G_{3,5} + G_{3,6}$, $\overrightarrow {G} = \left[ { - G_{3,1} , - G_{3,2} ,G(3,3), - G_{3,4} , - G_{3,5} , - G_{3,6} } \right]$ .c Illustration of a source node. Each numbered circle in the figure represents a node and each arrow represents a directed edge. The figure illustrates a graph segment with a source node (marked as s) and two vessels (GS1 and GS2). The resulted G matrix row vector for node s = 0 is: $\overrightarrow {G} = \left( { - 1, - G_{s,1} , - G_{s,2} } \right) = \left( { - 1, - G_{0,1} , - G_{0,2} } \right)$ , which results in the overall Nodal equation: $q_{s} = - (G_{s,1} + G_{s,2} )P_{0}$. d Illustration of a supply node. Each numbered circle or square in the figure represents a node and each arrow represents a directed edge. Dashed edges correspond to supply areas that are governed by auto-regulation. The figure illustrates two supply node—numbered 1 and 2. Sink represents the analogy to current sink. The resulted G matrix row vector for the supply node 1 is: $\overrightarrow {G} = - ( - G_{s,1} ,G_{s,1} )$, which results in the overall nodal equation: $q_{1} = Q_{s,1} = - G_{s,1} \cdot CPP$, where CPP stands for cerebral perfusion pressure

The mass flow conservation can be then applied using this equation for each node and the resulting set of equations can be solved to yield the pressures. In the presented application, it is applied via a modified nodal analysis (MNA) [21] that incorporates constraints on the system to be driven by a system pressure and ensures constant blood supply to specific regions.

Our network consists of three types of nodes and two types of edges. The primary type of node is one that connects different vessels with each other (i.e., junctions between vessels). The second type of node is a source that provides the system with blood. The third type of node is a supply area, whose resistance depends on the incident pressure and the auto-regulation process described in the next section (see Fig. 5). The types of edges in our network represent normal arterial blood vessels and vessels that connect the supply areas to the sink node. The resistance of the latter is determined via the auto-regulation function given in the next section.

Modeling auto-regulation

The simulated vessel network encompasses the circle of Willis and the larger arterial segments of 1st and 2nd order of the three major brain arteries: anterior-, medial- and posterior cerebral artery (ACA, MCA and PCA). Those are the vasculature segments that are (a) accessible for interventions such as surgery or thrombectomy and (b) their anatomical architecture can be derived from medical imaging. The vascular downstream regions after the above segments including the 3rd order segments represent a network of small arteries, finer small arterioles and the capillary bed that can change their radii in order to decrease or increase the blood supply. This process is called auto-regulation and ensures that the blood supply to the brain remains largely constant within certain limits. Most of the vascular network’s resistance originates from these supply areas. Autoregulation was implemented into the simulation framework according to the following equation based on literature values [36]:

$$ R_{{{\text{auto}}}} (P_{inp} ,P_{{{\text{out}}}} ) = \left\{ {\begin{array}{*{20}l} {\frac{{P_{{{\text{inp}}}} - P_{{{\text{out}}}} }}{{Q_{{{\text{supply}}}} }}} \hfill & {{\text{if}}\;150{\text{mm Hg}} > P_{{{\text{inp}}}} > 50{\text{mm Hg}}} \hfill \\ {R_{{{\text{min}}}} = \frac{{50{\text{mm Hg}}}}{{Q_{{{\text{supply}}}} }}} \hfill & {{\text{if}}\;P_{{{\text{inp}}}} < 50{\text{mm Hg}}} \hfill \\ {R_{{{\text{max}}}} = 2.25*R_{{{\text{min}}}} } \hfill & {{\text{if}}\;P_{{{\text{inp}}}} > 150{\text{mm Hg}}} \hfill \\ \end{array} } \right., $$

(6)

where $R\_auto$ is the autoregulated resistance of the vessel, $P\_inp$ and $P\_out$ are the input and output blood pressures of the vessel and $Q\_supply$ is the blood flow supply to the vessel.

In more detail, the blood flow into the supply areas is provided by the outgoing segments (A2, P2 and M1 or M2) of the circle of Willis, see Fig. 5. The behavior of the supply areas is modeled according to the autoregulation as described in Eq. 6. This means that the peripheral resistance of the supply area adjusts itself such that the blood flow is kept constant for a given pressure gradient $\Delta P$.

Boundary conditions

The model requires a systemic mean arterial pressure (MAP) that drives the flow through the network such that the blood supply to the various supply areas of the brain is constant, while the blood flows and pressures are adjusted accordingly. However, the auto-regulation model is limited to the region between 50 and 150 mmHg. Blood pressure below or above these boundaries indicates a pathological state in which the body cannot maintain the necessary blood pressure to ensure constant blood supply and results in hypoperfusion of the brain tissue. Under these conditions, the simulation indicates that the respective supply areas are not sufficiently supplied with blood. Clinically speaking, these areas are vulnerable to ischemia. In addition to the system pressure a blood flow rate per supply area was provided as P_in being the system pressure, while Q_s,i for i = {1,…,N} being the blood supply demand of supply area i as physical quantity.

Algorithmic derivation of the blood supply

The blood supply to the specified brain areas is derived using an adjusted modified nodal analysis (AMNA). The construction of the matrix equations per type of node and edge is detailed in the following. The described system is overdetermined by N equations, where N stands for the number of nodes. The last node is taken as the sink node, with a pressure value of 0. Therefore, the system is described by N−1 mass flow equations. These mass flow equations can be written in terms of a matrix:

$$ G \cdot \vec{P} = \vec{Q}, $$

(7)

where the matrix $G$ represents the network structure and consists of the conductivity (e.g., inverse resistances) of the vessels. $\overrightarrow {P}$ is the unknown blood pressure vector and $\overrightarrow {Q}$ is the blood flow vector. The AMNA algorithm yields a reduced size of the solution system by removal of known values from the solution vector $\overrightarrow {Q}$. The following section details the determination of the conductivity values and blood flow to construct the matrix equations from a vessel graph as depicted in Fig. 5 and followed by the AMNA implementation.

Junction nodes

The matrix row vector for a junction node is defined as:

$$ G(i,j) = \left\{ {\begin{array}{*{20}l} { - G_{i,j} , \quad i \ne j} \hfill \\ {\sum\limits_{j} {G_{i,j} } , \quad i = j} \hfill \\ \end{array} } \right., $$

(8)

where the diagonal element is the sum of all G values for all connected edges and the non-diagonal elements are set to be the negative conductivity value of the corresponding edge—or otherwise 0. An exemplary derivation of a junction node vector is described under Fig. 6b. The corresponding element of the $\overrightarrow {Q}$ is 0. This represents a flow equation as depicted by Eq. 5.

Source nodes

The matrix row vector for a source node is defined as:

$$ G(i,j) = \left\{ {\begin{array}{*{20}l} { - G_{i,j} ,\quad i \ne j} \hfill \\ { - 1 ,\quad i = j} \hfill \\ \end{array} } \right.. $$

(9)

The diagonal element is $G\left( {i,j} \right) \, = - 1$ and all other matrix elements are the negative conductivity values of the (connected) corresponding edges and otherwise 0, similarly to junction nodes. An exemplary derivation of the vector is described in Fig. 6c. The flow into the source node is determined via the supply areas and according to the mass flow can be written as:

$$ q_{s} = - P_{0} \cdot \sum {G_{s,i} } , \quad for \; Node \; i \; connected \; to \; source \; s. $$

(10)

For the nodes that are directly connected with the source node, the q vector element is then:

$$ q_{connected} = G_{s,connected} \cdot P_{0} . $$

(11)

Figure 6c depicts a graph segment with a source node and two connected vessels.

Supply nodes

To model supply areas, the auto-regulation equation is applied (Eq. 6). The matrix row vector for a supply node is defined as:

$$ G(i,j) = \left\{ {\begin{array}{*{20}l} { - G_{i,j} , \quad i \ne j} \hfill \\ {\sum\limits_{j} {G_{i,j} } , \quad i = j,j \ne sink} \hfill \\ \end{array} } \right.. $$

(12)

Similarly to junction nodes, the diagonal element $G\left( {i,j} \right)$ is the conductivity sum of incident edges, however in case the edge is connected to a sink it does not contribute to the sum. An exemplary derivation of the vector is described in Fig. 6d. The q-vector values are:

$q_{i} = - Q_{s,i}$ for the supply node i connected to source s.

And $q_{\sin k} = - Q_{s,i}$ for the sink connected to supply node i.

Figure 6d depicts a case with two supply areas, denoted by the thick dashed lines that are connected to the same sink.

Summary of nodal analysis construction

To summarize, the q-vector has a contribution for a supply area or a source node. If an edge connects a source with a supply area, the corresponding q-vectors will have contributions from both the source and the supply area. Table 3 summarizes the terms for the q-vector elements in dependence of the node type:

$$ G_{i,j} = \left\{ {\begin{array}{*{20}l} {R_{ext}^{ - 1} {\mkern 1mu} ,\;vessel} \hfill \\ {0{\mkern 1mu} ,\;\sup ply} \hfill \\ {R_{auto}^{ - 1} {\mkern 1mu} ,\;\sin k} \hfill \\ \end{array} } \right.. $$

(13)

Table 3 A: Summary of q vector element derivation in the modified nodal analysis (MNA)

Here vessel refers to a connection of junction node to junction node, supply refers to a connection of junction node to a supply node and sink refers to a connection of a supply node to sink. $R_{{{\text{ext}}}}$ and $R_{{{\text{auto}}}}$ are described in formulas 4 and 6, respectively.

Adjusted modified nodal analysis (AMNA)

The AMNA allows to reduce the size of the solution system by removing known values from the solution vector $\overrightarrow {Q}$ in consecutive 3 steps. Once the MNA matrix equations are determined, all the values of the system matrix that are known given system pressures are drawn to the right-side solution vector $\overrightarrow {Q}$. This pertains the values associated with the input pressures of the source nodes, i.e., MAP. As a result of this transition the equations associated with the source nodes become zero. In a second step, these redundant rows and columns are then deleted. Finally, the G matrix columns are swapped to yield a diagonal matrix, with the solution vector adjusted accordingly. This process simplifies the solution derivation and therefore accelerates its application.

Sensitivity analysis

In order to predict the effects of changes of a certain variable on the system, in this case the radius of a blood vessel, we performed a sensitivity analysis. Sensitivity analysis quantifies this effect by estimating the partial derivative of a system variable such as the blood pressure in this case, with respect to the radius of a given vessel. This is achieved by application of Newton–Raphson method using Taylor series expansion [37] as follows where—as previously defined—Q is the blood pressure and l is the length of the artery. The formula presented shows the sensitivity analysis—as an example—for the internal carotid artery denoted as int.car.I.:

$$ \frac{{\partial Q_{A1} }}{\partial rInt.Car.I} \approx \left[ {Q_{A1} (r_{A1} ,l_{A1} ,...,r_{Int.Car.I} ,l_{Int.Car.I} ,p_{in} ) \, - Q_{A1} (r_{A1} ,l_{A1} ,...,r_{Int.Car.I} + \Delta r_{Int.Car.I} ,l_{Int.Car.I} ,p_{in} )]/\Delta r_{Int.Car.I} } \right]. $$

(14)

Similar Newton steps can be formulated for all relevant system variables and allow to estimate how the system reacts to changes of certain system variables.

Simulation interface

The simulation interface was implemented as a java application with an integrated graphical user interface (GUI) under the loose-coupling paradigm to ensure that components can be exchanged easily. The key element of the simulation is a 2D projection as a representation of the simulated vessel tree. The interfacé́s main components are the simulation view including areas at risk and pressure view, and the possibility to change blood pressure as a boundary condition.

Perfusion imaging processing

DSC imaging was processed using the pgui software (Version 1.0, Center for functional neuroimaging, Aarhus University). Four arterial input functions were placed contralateral to the stenosis/occlusion in the M2 vessel area and visually assessed for optimal shape [16]. Deconvolution was performed according to the parametric method introduced by Mouridsen et al. [38]. Non-deconvolved time-to-peak, and deconvolved cerebral blood flow (CBF), time-to-maximum (Tmax) and mean-transit-time(MTT) maps were created and assessed in this study.

Comparison of simulation and perfusion imaging results

In 67 patients, we compared simulation and DSC perfusion imaging results. For this purpose, we defined vulnerability for both modalities and determined the detection rate of the simulation to identify vulnerabilities in DSC imaging by sensitivity and specificity.

Vulnerability of brain tissue to ischemia was defined in MR-imaging as a visually rated TTP and/or Tmax increase and/or CBF decrease. For the simulation, vulnerability was defined as a perfusion pressure below 50 mmHg, at normal blood pressure or at a mean-arterial-pressure of 70. We assessed the sensitivity/specificity of the simulation results to detect vulnerability as defined by the visual DSC-analysis. Results were recorded for anterior/middle/posterior cerebral artery regions (left/right).