Patient-Specific Modeling Could Predict Occurrence of Pediatric Stroke

Moyamoya disease (MMD) is a progressive steno-occlusive cerebrovascular disease leading to recurrent stroke. There is a lack of reliable biomarkers to identify unilateral stroke MMD patients who are likely to progress to bilateral disease and experience subsequent contralateral stroke(s). We hypothesized that local hemodynamics are predictive of future stroke and set out to noninvasively assess this stroke risk in pediatric MMD patients. MR and X-ray angiography imaging were utilized to reconstruct patient-specific models of the circle of Willis of six pediatric MMD patients who had previous strokes, along with a control subject. Blood flow simulations were performed by using a Navier–Stokes solver within an isogeometric analysis framework. Vascular regions with a wall shear rate (WSR) above the coagulation limit (>5,000 s⁻¹) were identified to have a higher probability of thrombus formation, potentially leading to ischemic stroke(s). Two metrics, namely, “critical WSR coverage” and “WSR score,” were derived to assess contralateral stroke risk and compared with clinical follow-up data. In two patients that suffered a contralateral stroke within 2 months of the primary stroke, critical WSR coverages exceeding 50% of vessel surface and WSR scores greater than 6× the control were present in multiple contralateral vessels. These metrics were not as clearly indicative of stroke in two additional patients with 3–5 year gaps between primary and contralateral strokes. However, a longitudinal study of one of these two cases, where a subsequent timepoint was analyzed, suggested disease stabilization on the primary stroke side and an elevated contralateral stroke risk, which was confirmed by patient outcome data. This indicates that post-stroke follow-up at regular intervals might be warranted for secondary stroke prevention. The findings of this study suggest that WSR-based metrics could be predictive of future stroke risk after an initial stroke in pediatric MMD patients. In addition, better predictions may be possible by performing patient-specific hemodynamic analysis at multiple timepoints during patient follow-up to monitor changes in the WSR-based metrics.

Introduction

Moyamoya disease (MMD) is a progressive cerebrovascular disorder characterized by stenotic or occlusive lesions in the terminal internal carotid artery (ICA) and the proximal middle or anterior cerebral arteries (MCA and ACA), leading to recurrent transient ischemic attacks and stroke. A substantial number of patients presenting with unilateral lesions go on to form bilateral disease, and it is known to occur more commonly in pediatric patients (Hirotsune et al., 1997). Neurosurgical intervention on the symptomatic arterial lesion followed by close observation for bilateral involvement is usually recommended (Hayashi et al., 2010), as unilateral stroke is sometimes followed by contralateral stroke(s), even after receiving surgical treatment. In our retrospective review of data from 50 pediatric MMD patients, we found that 39.2% of the unilateral stroke patients, 82% of whom underwent ipsilateral dural inversion surgery, suffered a subsequent contralateral stroke (Hossain et al., 2021). There is a lack of reliable biomarkers to identify these at-risk patients and intervene before the subsequent stroke occurs. Progression of vasculopathy associated with MMD is highly variable from patient to patient, and patients with similar angiographic appearance can have varying degree of stroke risk (Zipfel et al., 2009). As a result, several studies that have relied on angiographic appearance to predict progression to bilateral disease, produced conflicting conclusions (Yeon et al., 2011).

While pathophysiology of MMD remains unclear, hemodynamic stress is known to play an important role in stroke (Derdeyn, 2001; Zipfel et al., 2009). A vascular wall shear rate (WSR) above the coagulation limit (>5,000 s⁻¹) signifies an increased probability of thrombus formation (Casa et al., 2015), leading to stroke (Kim et al., 2017). In a preliminary study using an Acta2^−/− mouse model that develops many of the features of MMD (Milewicz et al., 2010; Starosolski et al., 2015), we demonstrated that occlusion in one of the major arteries in the circle of Willis (CoW) increased hemodynamic burden in vessels contralateral to the occlusion, thus potentially raising contralateral stroke risk if critical WSR limits are exceeded (Hossain et al., 2021). If a similar behavior is realized in the MMD patients, it could have a profound implication for patient care. We therefore hypothesized that local WSR could be predictive of future stroke and that image-based patient-specific hemodynamic analysis, in conjunction with clinical observations, could be used to noninvasively assess this stroke risk in pediatric MMD patients.

In a recent study, Lee et al. hypothesized that endothelial shear stress could act as a trigger for alteration of vascular remodeling and correlated contralateral progression of unilateral MMD to spatial variability of endothelial shear stress around contralateral terminal ICA (Lee et al., 2020). This wall shear stress was indirectly estimated by using signal intensity gradient (SIG) in time-of-flight (TOF) sequences from brain MRI, but validation with computational fluid dynamics (CFD) was not performed. In the present work, we conduct physiologically realistic CFD analysis of the full CoW and directly quantify local WSR distribution for a pilot cohort (n = 6) of pediatric MMD patients. The objective of this retrospective study is to identify susceptible vascular regions with WSR critically above 5,000 s⁻¹ that could evolve into severe stenosis or complete occlusion due to thrombus formation, leading to ischemic strokes (Casa et al., 2015). The in-silico results are compared with patients’ outcome data to investigate whether such subject-specific analysis of local WSR utilizing patient’s post-stroke MR TOF and X-ray angiography (XA) images is predictive of future contralateral stroke. The findings of this work could provide the necessary basis for a more comprehensive study involving a large cohort of patients that could inform follow-up strategies for unilateral patients at risk for contralateral stroke, and potentially guide timing of preventive interventions.

Materials and Methods

Vascular Model Creation From Imaging Data

The generation of patient-specific CoW models from MR TOF images and their adjustment with XA data are outlined below.

MR TOF Segmentation

Patient MR TOF images (Supplementary Figure S1A) are segmented in 3D Slicer using a combination of intensity thresholding and manually painting vessels in each slice. The final segmented label map is exported as a triangulated 3D surface mesh and is considered the initial model (Supplementary Figure S1B).

Diameter Adjustment Using XA Images

3D MR TOF images can have poor signal where blood flow is slow or complex, causing vessel caliber to be inaccurately represented (Hoogeveen et al., 1998). Inadequate image resolution can cause areas of MMD-related severe narrowing to be poorly resolved, resulting in inaccurate evaluation of local hemodynamics (Hossain et al., 2021). To improve accuracy of vessel geometry (e.g., vessel diameter), the MR-derived CoW model (initial model) is adjusted based on 2D imaging data obtained through XA (Supplementary Figure S1C), which is considered the gold standard for evaluating vessel narrowing.

Details of this method can be found elsewhere (Horn et al., 2021b). Briefly, the initial model is first aligned to the XA imaging view to enable comparison of the initial CoW model vessel geometry with that seen on 2D XA images. For each XA image, an imaging axis is defined in a virtual space by using XA imaging metadata, including “distance source to patient,” “distance source to detector,” “positioner primary angle,” and “positioner secondary angle.” Centerlines are extracted from the initial model and projected along the imaging axis into the virtual detector plane. Using custom algorithms implemented in MATLAB (version 2019b, The MathWorks Inc., Natick MA), alignment of the centerlines with the XA image is ensured. This enables direct comparison between the clinical images and the initial CoW model.

Next, algorithmic adjustments to the initial CoW model are made based on each XA image using 3D computer-aided design (CAD) software Rhino and Grasshopper as follows. The vessel diameters of the CoW model are measured at various locations along each vessel and projected into the virtual detector plane of XA image. The vessel edges are manually drawn on the XA images and target vessel diameters are obtained at those locations. The model vessel diameters are then compared to the target vessel diameters to compute scale factors, which are used to locally adjust model vessel diameter radially about the vessel centerline.

Finally, to achieve a better agreement with XA images, iterative manual adjustments are made to the 3D model geometry. In this procedure, virtual angiographies of the 3D adjusted model are projected into the virtual detector plane of each XA image using the alignments obtained earlier. The resulting 2D projections of the model are overlaid into the XA images for comparison. In regions, where the 2D projection does not match the vessel geometry as seen on the XA image, the model vessel geometry (i.e., vessel diameter and vessel curvature) is manually adjusted using MeshMixer (version 3.5, Autodesk Inc., San Rafael, CA). This adjustment-projection-comparison loop is iteratively repeated until the 2D projections of the model and the corresponding XA images agree, as confirmed by a neuroradiologist. The final adjusted model (Supplementary Figure S1D) is then exported as a triangulated surface mesh.

Reconstruction of Solid NURBS Mesh

The adjusted surface mesh is used to generate a volumetric non-uniform rational B-spline (NURBS) reconstruction of the CoW to be used in analysis. The methodology, detailed in a previous report (Urick et al., 2019), uses an in-house template-based vascular modeling software that leverages the geometric modeling kernel of the CAD package Rhino for its robustness, interoperability, accuracy, and speed. First, centerlines are extracted from the triangulated surface mesh with a skeletonization algorithm. The centerlines are sampled at regular intervals, and at each sample point, a perpendicular frame is constructed. At branch points, where multiple vessel centerlines meet, the frames are folded to build a conforming parameterization. A mesh-frame intersection curve is computed at each frame and interpolated with a NURBS curve. The NURBS curves are then lofted together along the centerlines to form a NURBS reconstruction of the triangulated surface mesh geometry that is parameterized in the circumferential and axial directions. Finally, a radial direction is added to the parameterization by extruding the surface control points to the centerline, yielding a hexahedral NURBS solid model (Supplementary Figure S1E) that is used for analysis. Mesh sizes and refinement strategies were chosen based on mesh independence studies (Supplementary Figure S2; Hossain et al., 2021).

Governing Equations and Solution Strategies

Blood flow was assumed to be governed by the unsteady Navier–Stokes equations with a time-dependent pulsatile inflow boundary condition (Wahlin et al., 2013) imposed at the three inlets: left ICA (LICA), right ICA (RICA), and basilar artery (BA), a no-slip boundary condition was prescribed at the rigid vascular wall, and a traction-free outflow boundary condition implemented at the branch outlets (see Figure 1 for the simulation setup). Blood was modeled as an incompressible Newtonian fluid with a density ρ of 1,060 kg/m³ and a dynamic viscosity μ of 0.0035 Pa∙s. A Navier–Stokes solver within an isogeometric analysis framework (Hsu et al., 2011; Hossain et al., 2012, 2013) was used to solve the system of equations by applying the solution strategies and numerical procedures described in a previous work (Hossain et al., 2021). The WSR was computed by using the following equation of wall shear stress vector $\tau =\left(\left(,\sigma ·n,\right)\right)-\left(\left(\left(,\right(,\sigma ·n,)\right)·n,)\right)n$ and the relation $\mathrm{WSR}=\frac{1}{\mu }\left(\left|,\tau ,\right|\right)$ , where $\left(\left(,\sigma ·n,\right)\right)$ is the traction vector, $\sigma$is the stress tensor, and $n$is the unit normal.

FIGURE 1

Figure 1. Simulation setup. (A) The hexahedral NURBS mesh of a human CoW with a bifurcation (B) is shown. A cross-section of the computational mesh (C) is shown to highlight the boundary layer refinement (D) we implemented. The governing equations subjected to boundary conditions (B.C.s) are also presented, where u represents velocity, p is the pressure, f is the external force, μ is dynamic viscosity, ρ is density, t is time, and n is the unit normal. (E) A pulsatile inflow condition is prescribed at the three inlets. Left and right ICAs (red waveform) and BA (blue waveform). A no-slip B.C. is imposed on the vessel wall, and a traction-free B.C. is set at all the outlets.

Results

Data from 50 pediatric MMD patients were retrospectively reviewed (see Supplementary Figure S3; Hossain et al., 2021) and six cases were selected for CFD analysis (Table 1) based on the following criteria: (1) patient age is less than 18 years, (2) no other cerebral or cerebrovascular disease is indicated, (3) history of at least one stroke, (4) availability of post-stroke MR TOF and XA imaging taken within 6 months of each other, and (5) there is communication between posterior and anterior circulations in the CoW. For each MMD patient, blood flow velocity fields were simulated in the CoW models reconstructed from imaging data and WSR values were quantified. The spatial distribution of WSRs at peak systole is shown in Figure 2, along with timelines summarizing the chronology of stroke events, surgeries, and the collection of imaging data used to generate each model. The results are compared with those in a pediatric control subject (see Supplementary Figure S4) presented in a previous report (Hossain et al., 2021).

TABLE 1

Table 1. Summary of the MMD patient cohort (n = 6).

FIGURE 2

Figure 2. Distribution of WSR above the coagulation limit (>5,000 s⁻¹) at peak systole (top) along with corresponding patient timeline (below) for the patients we studied. (A) MMD1, (B) MMD2, (C) MMD3, (D) MMD4, (E) MMD5, and (F) MMD6. Red circles indicate stroke occurrence, green circles indicate surgery, yellow ovals denote model creation from XA and MR TOF imaging data, and red circle octagon is the last follow-up imaging available.

WSR Distribution and Patient Outcome

Patient MMD1 presented with right-sided stroke and was treated with right dural inversion surgery. Post-operative images were used for CoW model creation and CFD analysis (Figure 2A). At peak systole, simulation results predict extreme WSR levels in areas of vessel narrowing in right MCA (RMCA), right ACA (RACA), and RICA, (maximum 64,176 s⁻¹, 61,294 s⁻¹, and 44,617 s⁻¹, respectively), coinciding with the location of the primary stroke. Concurrently, critical WSRs were seen in the terminal LICA (maximum 63,094 s⁻¹). The patient experienced a contralateral stroke 2 months later.

Patient MMD2 initially presented with a right-sided stroke, followed by a left-sided stroke approximately 2 years later. A right dural inversion surgery was performed shortly after the second stroke. Post-operative images were used for model creation. At peak systole, WSR levels in the RACA and RMCA (Figure 2B) reached 33,728 s⁻¹ and 40,296 s⁻¹, respectively, coinciding with areas with stenosis. Severe narrowing was also observed in the left ACA (LACA) with a maximum WSR of 15,796 s⁻¹. While not visibly impacted by MMD-related narrowing, the left MCA (LMCA) had elevated WSRs (maximum 28,370 s⁻¹). A left dural inversion surgery was performed within 1 year of the left-sided stroke, and the disease stabilized during 7 year follow-up, at the end of which, evidence of bilaterally diffuse ischemic injury was noted. It is unclear from the patient history if this is related to the previously mentioned bilateral stroke events.

Patient MMD3 suffered two left-sided strokes within the same year. Post-stroke images were used to create a CoW model. Extreme WSRs at peak systole were observed (Figure 2C) in the LACA and LMCA (maximum 27,256 s⁻¹ and 45,066 s⁻¹, respectively), coinciding with the primary stroke location. Contralateral WSRs were considerably lower, peaking around 12,761 s⁻¹ in the RMCA. Left side and right-side dural inversion surgeries were performed within one and 2 years of the primary stroke, respectively. A right-sided stroke occurred approximately 3 years later.

Patient MMD4 presented with a right-sided stroke. Post-stroke images were utilized for model creation. Super critical WSR vales at peak systole were predicted (Figure 2D) in the severely narrowed RACA (maximum 72,588 s⁻¹) and in the RMCA (maximum 45,831 s⁻¹) where an ischemic event was evident. Concurrently, elevated WSR levels were seen in the LMCA (21,203 s⁻¹) and, despite showing vessel narrowing, the LACA showed no WSRs at peak systole above the coagulation threshold. Within a year of the primary stroke, patient underwent a right-side surgery, which was closely followed by two contralateral strokes. After a subsequent surgery on the left side, the disease stabilized.

Patient MMD5 suffered a left-sided stroke and within a year underwent left-side surgery. Post-operative images were utilized to create a CoW model. Elevated peak systolic WSR values (maximum 30,974 s⁻¹) were predicted (Figure 2E) in the severely narrowed region of the LACA and the LMCA (maximum 19,305 s⁻¹). Contralaterally, the RMCA and the RACA, the latter of which showed significant stenosis, saw peak WSR values of 12,015 s⁻¹ and 11,829 s⁻¹, respectively. Within 3 years of the primary stroke, the patient suffered bilateral strokes.

Patient MMD6 presented with a left-sided stroke. Model created using post-stroke images shows that peak WSRs reached 22,140 s⁻¹ in the LMCA, where there is a mild narrowing, and 18,060 s⁻¹ in the terminal LICA (Figure 2F). Contralaterally, the RMCA and terminal RICA had similarly elevated peak WSRs (17,865 s⁻¹ and 20,980 s⁻¹, respectively). Within 1 year of the primary stroke, bilateral dural inversion surgery was performed. The patient remained stable throughout the next 10 years of follow-up.

WSR-Based Metrics

To further analyze the impact of elevated WSR values on stroke risk, we generate WSR distribution curves as described in Figures 3A,B for the right and left MCAs, ACAs, and terminal ICAs, where the terminal ICA is defined as the ICA segment between the ophthalmic artery and the ICA bifurcation. The WSR distribution curve for a vessel shows how much of the vessel’s surface area (SA) exceeds a given WSR threshold value. Figures 3C,D illustrate these curves for the left side vessels of the control subject and patient MMD3, respectively, indicating that vessel areas of higher WSR are more widespread for patient MMD3 compared to the control. From the WSR distribution curves, two metrics are derived to assess stroke risk.

FIGURE 3

Figure 3. (A) Construction of WSR curve. First, a vessel segment is extracted from the model and the total surface area is measured. A WSR threshold is set and the surface area that exceeds this threshold is quantified and then divided by the vessel’s total surface area to obtain percent vessel area exceeding WSR threshold, which is then calculated for multiple threshold values ranging from 0 to 100,000 s⁻¹ with an increment as low as 100 s⁻¹ to generate the WSR curve (B). The WSR curves for the left and right ICAs, MCAs, and ACAs are shown for (C) the control subject and (D) MMD3.

First, we consider the “critical WSR coverage” (CWSRC) in each vessel representing the percent of each vessel’s surface area with critical WSR values—i.e., those exceeding the coagulation threshold of 5,000 s⁻¹ (Figure 4). A graphical explanation of how this metric is derived from the WSR distribution curve is shown in Supplementary Figure S5. In the control, three vessels (LMCA, LACA, and RACA) had little to no surface area with critical WSRs (0.9, 0.0, and 0.0%, respectively). The remaining vessels under consideration have mild to moderate CWSRC (LICA: 15.8%, RMCA: 29.0%, and RICA: 35.9%). In contrast, many of the vessels in the six MMD patients exhibit greater CWSRC than those in the control. Of note, the RMCA of MMD1 saw critical WSRs covering 88.3% of its surface area, and MMD4 saw high CWSRC in the LMCA (77.1%) and the RMCA (96.5%). As an aggregate form of this metric, the three vessels on either side are combined to form two regions, left and right, in which the CWSRC is assessed (Figure 5). For the control, the aggregate CWSRC values were 6.3 and 20.5% for the left and right sides, respectively. In the MMD patients, these aggregates vary from 22 to 52.2% on the left side and 13.7 to 63.7% on the right side, and each patient has at least one side with an aggregate above 20.5%, the greater of the two aggregates in the control.

FIGURE 4

Figure 4. Critical WSR coverage, expressed as a percentage of each vessel’s surface area with WSR >5,000 s⁻¹, in the left and right MCAs, ACAs, and ICAs for the control subject and the six MMD patients we studied. The label (Stroke) indicates the side of primary stroke.

FIGURE 5

Figure 5. Aggregate critical WSR coverages for the left and right sides of the anterior circulation (MCA, ACA, and ICA) for the control subject and the six patients we studied as: MMD1, MMD2, MMD3 and MMD3b, MMD4, MMD5, and MMD6. The percentage is calculated as the sum of the surface area exceeding WSR >5,000 s⁻¹ in the three vessels divided by the sum of the surface areas of the three vessels. The label (Stroke) indicates the side of primary stroke event.

The second metric we examine is a “WSR score” for each vessel calculated by integrating the WSR distribution curve over the domain (5,000, ∞) above the coagulation limit (see Supplementary Figure S6). This parameter gives an indication of the extent and magnitude of the higher WSRs. For example, two vessels with the same CWSRC can have different WSR scores if one of them has higher WSR values occupying the same vessel surface area. WSR score is relatively small or zero for the vessels in the control while many of the vessels in the MMD patients have significantly higher WSR scores with each having at least two vessels with a WSR score at least 4× the maximum seen in the control, and all but one (MMD6) with at least one vessel with a WSR score greater than 8× the maximum control WSR score (Figure 6). The MMD patients have left- and right-side aggregate WSR scores ranging from 1.5× to 37.5× the control’s largest aggregate WSR score (Figure 7).

FIGURE 6

Figure 6. WSR score in the left and right ICAs, MCAs and ACAs for the control subject and the six MMD patients we studied. The label (Stroke) indicates the side of primary stroke.

FIGURE 7

Figure 7. Aggregate WSR score. Area under the WSR distribution curves from Figure 4 aggregated by side for the control subject and the six patients we studied as: MMD1, MMD2, MMD3 and MMD3b, MMD4, MMD5, and MMD6. For each patient, the values shown in Figure 5 are added. The label (Stroke) indicates the side of primary stroke event.

A Longitudinal Study

Figure 8 summarizes a longitudinal study conducted for patient MMD3 where a second CoW model (MMD3b) was generated using imaging data collected soon after the second (right-side) dural inversion surgery to assess how disease progression impacts the WSR-based metrics for predicting future stroke risk. The most notable change is that the LACA has become fully occluded. Additionally, corresponding with approximately 2 years of patient growth, the vessels of the second model generally have larger diameters. As a result, the WSR values predicted at peak systole are less extreme for MMD3b (global maximum WSR of approximately 17,000 s⁻¹) compared to MMD3 (45,000 s⁻¹). On the left side, marked reductions in both metrics were seen compared to the initial model MMD3 (Figures 8H,I). CWSRC decreased in the LMCA (down to 12.9% from 52.2%) and in the terminal LICA (29.1% versus 58.5%), and the left side aggregate CWSRC decreased from 52.1 to 22.0%. The WSR scores in the LMCA and LICA also decreased by 87 and 71%, respectively, and the aggregate WSR score on the left decreased by approximately 88%. Contralaterally, a larger surface area of the RMCA and RACA exceeded the coagulation threshold (from 38.7 to 46.8% and from 13.1 to 28.9%, respectively), and the aggregated CWSRC on the right increased from 17.6 to 18.7%. The WSR score in the RMCA and RACA was approximately 3× and 17× that of the first model, respectively, and the right-side aggregate WSR score was approximately 3.7× that of the first model.

FIGURE 8

Figure 8. A representative case. Patient presented with multiple transient ischemia attacks on the left side. Initial brain MR TOF (A) and XA (B) confirmed occlusion of left M1 and left A1, and mild to moderate narrowing of right A1 and right M1 segments. The patient underwent unilateral left dural inversion surgery. The follow-up MR TOF (C) and XA (D) imaging within 2 years from the initial diagnosis show evidence of moderate contralateral progression of disease and right dural inversion surgery. (E) In the follow-up transverse (axial) soft tissue MR image (SENSE sequence) at approximately 5 years from the initial diagnosis, evidence of right-sided stroke within perirolandic region involving right MCA territories. (F) WSR distribution pattern in the patient’s CoW model MMD3 created from post-operative images shown in A and B, and (G) WSR distribution pattern for CoW model MMD3b created from post-operative images shown in C and D. Plots reporting the critical WSR coverage (H) and WSR score (I) of each vessel for the second CoW model (MMD3b) with comparison to the results from MMD3 (gray bars).

Discussion

We performed blood flow analysis of patient-specific models of the CoW in pediatric MMD patients and showed that local hemodynamics are markedly different than those present in a healthy subject (Hossain et al., 2021). These differences, which are largely dependent on vessel morphology altered by MMD progression, may be indicative of future contralateral stroke risk in patients presenting with unilateral stroke. We hypothesized that a critically high WSR above the coagulation limit of 5,000 s⁻¹ could cause thrombus formation, potentially leading to an ischemic event (Kim et al., 2017) and a careful analysis of local WSR could be predictive of future stroke. As such, we focused on the spatial distribution of WSR of the major vessels of the anterior circulation: ICA, MCA, and ACA (Kelly et al., 2006) and considered two WSR-derived metrics. First, we computed CWSRC, the percentage of the vessel’s surface area with critical WSR values exceeding the coagulation threshold. Second, we calculated a WSR score that gives a sense of the magnitude and extent of the critically high WSR values present in a vessel.

Among the six patients, we selected for analysis, all of which had at least one unilateral stroke event prior to the collection of the images used for model creation, three (MMD1, MMD3, and MMD4) suffered contralateral stroke after being surgically treated on the affected side, one (MMD5) had bilateral stroke after being surgically treated on the affected side, one (MMD2) experienced bilateral stroke prior to the imaging used for model creation, and one (MMD6) had no additional strokes within the 10 years of follow-up after being treated with bilateral surgery. In the patient with a history of bilateral stroke (MMD2), both right and left MCAs have more than 50% CWSRC and WSR scores >9× the max WSR score of the vessels in the control. The metrics are similarly high in the other patients in the vessels ipsilateral to the primary stroke. One or more vessels with >50% surface area with critical WSRs and WSR scores >7× the maximum control WSR score are seen in all these cases, which is consistent with the stroke events the patients experienced.

In two patients (MMD1 and MMD4), contralateral strokes occurred shortly after the primary stroke and being surgically treated on the affected side. High metric values (more than 50% CWSRC and WSR scores greater than 6× the max control WSR score) were seen in at least two of the contralateral vessels, potentially indicative of the impending contralateral stroke events. As further evidence of contralateral stroke risk for MMD1, the aggregated CWSRC on the contralateral side (52.2%) is greater than that on the primary stroke side (41.4%), and the aggregate WSR score, which is approximately 22× the maximum aggregated WSR score in the control, is also 7.5% higher on the contralateral side. This may indicate that at the timepoint represented by this CoW model, the disease state on the left side may have become more severe than what caused the primary right-sided stroke. On the other hand, for MMD4, the aggregate CWSRC on the contralateral side is comparable to the ipsilateral side (50.6% vs. 63.7%, respectively), but the contralateral aggregate WSR score is only 24% of the ipsilateral side. However, this value is still approximately 9× that seen in the control, indicating that increased contralateral stroke risk may exist.

For the other two patients (MMD3 and MMD5) that experienced contralateral strokes, these metrics do not strongly suggest future contralateral stroke. There were no contralateral vessels with >50% CWSRC or with WSR scores >3× the max control WSR score. The contralateral aggregate CWSRCs are much less when compared to the stroke side for MMD3 (17.6% vs. 52.1%) and MMD5 (13.7% vs. 29.1%) and the aggregate WSR scores in vessels contralateral to the primary stroke side are 10 and 24% of the ipsilateral aggregates for MMD3 and MMD5, respectively. While none of the WSR metrics strongly indicate future stroke risk, in both cases, there is a multi-year gap between the imaging used for model creation and the subsequent contralateral stroke(s). This suggests that perhaps, we can reasonably predict future strokes in the near term and that additional modeling, at additional timepoints during follow-up, may be warranted for patients without conclusively high or comparable (to ipsilateral) metric values.

To assess how local hemodynamics changes over time as the disease progresses, we examined imaging data at a subsequent timepoint for patient MMD3, approximately 2 years after the primary stroke. Comparing the two modeling timepoints (i.e., MMD3 and MMD3b), both CWSRC and WSR score decrease for all three vessels ipsilateral to the primary stroke, possibly indicating that the disease is stabilizing on the left side. Conversely, marked increases in the two metrics are seen in the RMCA and RACA. The aggregate WSR score on the contralateral side is >3× the score on the left side, suggesting that the patient is at risk of a contralateral stroke. Patient MMD3 did indeed go on to have a right-sided stroke 3 years later. Overall, the results from this longitudinal study of patient MMD3 demonstrate the utility in regular patient monitoring and performing hemodynamic analysis using updated CoW models to assess changes in stroke risk over time.

Assessing multiple timepoints may have been particularly useful for patient MMD6. The CWSRC in the vessels contralateral to the primary stroke was comparable to the ipsilateral side (35.0% vs. 39.8%, respectively). While the same pattern holds for the WSR scores, aggregate WSR score is approximately 31% greater on the contralateral side. Although similar evidence in other patients was predictive of future contralateral stroke, no additional stoke events are noted in this patient’s history. Notably, this patient had a bilateral surgery shortly after the primary stroke, which is known to improve cerebral perfusion and help washout thromboemboli in MMD patients (Kim et al., 2017), thereby potentially preventing further stroke.

A major goal of this work is to assess clinical value to evaluating patient-specific WSR as a potential biomarker for assessing future ischemic stroke risk because cerebral angiographic appearance, including degree of vessel narrowing, has not proven to be consistently reliable (Yeon et al., 2011), as it has been reported that patients with similar angiographic appearance can have different stroke outcomes (i.e., stroke or no stroke in future; Zipfel et al., 2009). In the blood flow simulations of the pediatric MMD patients studied in this work, we observed that patients with similar vascular narrowing features can have very different hemodynamics throughout various portions of the circle of Willis, and vice versa. For example, in patient MMD1, who initially presented with a right-sided stroke, the maximum WSR seen on the ipsilateral side (64,176 s⁻¹ in RACA) coincided with the narrowest stenosis seen (approximate diameter 0.0402 cm) on the right side of the brain. The patient underwent right dural inversion surgery. At that time, the narrowest stenosis seen contralaterally was in the LICA (approximate diameter 0.141 cm) with a WSR of 40,379 s⁻¹. But the maximum WSR on the left side occurred at terminal LICA (63,094 s⁻¹) with a stenosis diameter of approximately 0.162 cm. If we were to rely only on the degree of stenosis, a left-sided imminent stroke would not have been suspected (which did occur 2 months later), as these vessels contralateral to the primary stroke were larger by a factor of 3 or 4. However, according to our analysis, although LICA diameter was significantly larger than RACA, local WSR was critically high (9× that seen in control) at the terminal LICA location and comparable to that seen on the stroke side (RACA). Hence, our recommendation would have been a bilateral surgery (instead of the unilateral surgery that this patient had) because our analysis suggested that the patient might be at risk for a contralateral stroke. We therefore believe that a quantitative assessment of local WSR distributions in all vessels of the CoW has added value in the prediction of future stroke occurrence because it is based on the hemodynamics-triggered mechanism of thrombus formation, which can lead to thrombotic or thromboembolic stroke. In a future larger study, it might be worth investigating if our WSR-based approach does indeed fare better than a simple measurement of stenosis.

In this work, we adopt a WSR of 5,000 s⁻¹ as a limit for coagulation, above which localized platelet accumulation can occur on a thrombogenic surface (Casa et al., 2015). However, in the MMD patients we studied, stroke events did not seem to occur immediately upon exceeding this coagulation limit. The lowest WSR that could be associated with a subsequent contralateral stroke was approximately 21,000 s⁻¹ (MMD4). In patient MMD1, even though the contralateral WSRs reached 12× the coagulation limit and 5× that seen in the control subject, contralateral stroke did not occur until 2 months later. The mechanisms of stroke in MMD are still unclear (Kim et al., 2017). Microfluidics models studying the effect of WSR on platelet rich thrombus formation suggest that pathologically high flow rates induce endothelial damage, platelet aggregation, and fibrin disposition (Tsai et al., 2012). In MMD, pathogenesis of stroke could be explained as follows. Upon thrombus formation, entrapped cells remain in the intimal layer even after the clot resolves. With repeated clot formation, intimal hyperplasia results and the lumen narrows. An acute stroke can result when a clot or embolus finally occludes the lumen. Further studies of MMD pathology are required to understand this phenomenon.

The study presented here analyzed a limited cohort selected from a single-hospital dataset, and thus, our findings are likely limited by selection bias and the inherent biases associated with a retrospective study. While the results presented in this work are not statistically conclusive due to the small sample size, we do see general trends that are interesting enough to warrant further investigation. Longitudinal studies with larger sample sizes are needed to confirm that the WSR-based metrics identified in this work are predictive of future stroke. This should include analysis of how timing of neurosurgical intervention impacts modeling results. In the present retrospective study, due to limited availability of imaging data, the models were created irrespective of when surgery occurred, if it occurred at all. Future prospective studies should also aim to measure flowrates in the ICAs and BA at time of imaging (e.g., including 4D MRI sequences to the MR TOF acquisition or utilizing carotid ultrasound to measure mean flow rates) so that patient-specific boundary conditions can be prescribed in the simulations.

In summary, we describe two possible WSR-based metrics for assessing stroke risk in pediatric MMD patients. In a subset of the patients evaluated, these metrics strongly suggest imminent contralateral stroke, which was confirmed by patient outcome data. Inconclusive results in the remaining patients highlight the difficulty in making predictions based on the hemodynamic environment at a single time point. Changes in the WSR-based metrics over time could provide a better assessment of future stroke risk. Given the limited number of patients studied herein, a prospective longitudinal study of a larger randomized patient cohort is warranted. If shown promise in the larger study, our patient-specific computational modeling approach could be used to noninvasively predict future stroke occurrences in pediatric MMD patients and help ensure that unilateral MMD patients are properly monitored and treated before strokes occur.

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.

Ethics Statement

The studies involving human participants were reviewed and approved by Institutional Review Board of Baylor College of Medicine and affiliate hospitals. Written informed consent from the participants’ legal guardian/next of kin was not required to participate in this study in accordance with the national legislation and the institutional requirements.

Author Contributions

JH performed vascular model segmentation and model adjustment, postprocessed, analyzed, and interpreted the data. MJ generated NURBS meshes from segmented models. JH and MJ prepared the Figures. ZS selected the clinical subjects for the study, performed preliminary image processing, oversaw model segmentation and model adjustment, and interpreted the results. AM reviewed 2D model projections, brain MRIs and angiograms to confirm model agreement with angiogram data, and performed Suzuki staging. DM advised on the clinical relevance of the work. AA interpreted results. JH and SH drafted the manuscript. SH conceptualized, planned and supervised the project, performed simulations, and interpreted results. All authors provided critical feedback, edited the manuscript, and approved the final submission.

Funding

This work was supported by NIH grant R03NS110442 awarded to SH. DM acknowledges support from AHA Merit Award and the Olivia Petrera-Cohen Research Fund.

Conflict of Interest

AA is a consultant to, and a founder and stockholder in, Alzeca Biosciences and a shareholder in Sensulin LLC. ZS is a stockholder in Alzeca Biosciences and a consultant for InContext.ia.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors gratefully acknowledge Karen Chen, MD at Texas Children’s Hospital for helping with Suzuki staging, Mario Dada at the University of Texas at Austin for assisting with image segmentation and Texas Advanced Computing Center (TACC) for providing high-performance computing resources. The content of the manuscript has previously appeared online in a preprint (Horn et al., 2021a).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.846404/full#supplementary-material

Abbreviations

MMD, Moyamoya Disease; CoW, Circle of Willis; WSR, Wall Shear Rate; ICA, Internal Carotid Artery; MCA, Middle Cerebral Artery; ACA, Anterior Cerebral Artery; BA, Basilar Artery; CAD, Computer Aided Design; NURBS, Non-Uniform Rational B-Splines; CWSRC, Critical Wall Shear Rate Coverage.

References

Casa, L. D., Deaton, D. H., and Ku, D. N. (2015). Role of high shear rate in thrombosis. J. Vasc. Surg. 61, 1068–1080. doi: 10.1016/j.jvs.2014.12.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Derdeyn, C. P. (2001). Hemodynamic impairment and stroke risk: prove it. AJNR Am. J. Neuroradiol. 22, 233–234.

PubMed Abstract | Google Scholar

Hayashi, K., Suyama, K., and Nagata, I. (2010). Clinical features of unilateral moyamoya disease. Neurol. Med. Chir. 50, 378–385. doi: 10.2176/nmc.50.378

CrossRef Full Text | Google Scholar

Hirotsune, N., Meguro, T., Kawada, S., Nakashima, H., and Ohmoto, T. (1997). Long-term follow-up study of patients with unilateral moyamoya disease. Clin. Neurol. Neurosurg. 99, S178–S181. doi: 10.1016/S0303-8467(97)00043-7

CrossRef Full Text | Google Scholar

Hoogeveen, R. M., Bakker, C. J. G., and Viergever, M. A. (1998). Limits to the accuracy of vessel diameter measurement in MR angiography. J. Magn. Reson. Imaging 8, 1228–1235. doi: 10.1002/jmri.1880080608

PubMed Abstract | CrossRef Full Text | Google Scholar

Horn, J. D., Johnson, M. J., Starosolski, Z. A., Meoded, A., Milewicz, D. M., Annapragada, A., et al. (2021a). Toward noninvasive assessment of stroke risk in pediatric cerebrovascular disease. medRxiv [Preprint]. Available at: http://medrxiv.org/content/early/2021/12/21/2021.12.17.21267944.abstract (Accessed December 21, 2021).

Google Scholar

Horn, J. D., Starosolski, Z. A., Johnson, M. J., Meoded, A., and Hossain, S.S. (2021b). A novel method for improving the accuracy of MR-derived patient-specific vascular models using X-ray angiography. bioRxiv [Preprint]. Available at: http://biorxiv.org/content/early/2021/12/23/2021.12.22.472309.abstract (Accessed December 23, 2021).

Google Scholar

Hossain, S. S., Hossainy, S. F. A., Bazilevs, Y., Calo, V. M., and Hughes, T. J. R. (2012). Mathematical modeling of coupled drug and drug-encapsulated nanoparticle transport in patient-specific coronary artery walls. Comput. Mech. 49, 213–242. doi: 10.1007/s00466-011-0633-2

CrossRef Full Text | Google Scholar

Hossain, S. S., Starosolski, Z., Sanders, T., Johnson, M. J., Wu, M. C. H., Hsu, M.-C., et al. (2021). Image-based patient-specific flow simulations are consistent with stroke in pediatric cerebrovascular disease. Biomech. Model. Mechanobiol. 20, 2071–2084. doi: 10.1007/s10237-021-01495-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hossain, S. S., Zhang, Y., Liang, X., Hussain, F., Ferrari, M., Hughes, T. J., et al. (2013). In silico vascular modeling for personalized nanoparticle delivery. Nanomedicine (Lond.) 8, 343–357. doi: 10.2217/nnm.12.124

PubMed Abstract | CrossRef Full Text | Google Scholar

Hsu, M.-C., Akkerman, I., and Bazilevs, Y. (2011). High-performance computing of wind turbine aerodynamics using isogeometric analysis. Comput. Fluids 49, 93–100. doi: 10.1016/j.compfluid.2011.05.002

CrossRef Full Text | Google Scholar

Kelly, M. E., Bell-Stephens, T. E., Marks, M. P., Do, H. M., and Steinberg, G. K. (2006). Progression of unilateral moyamoya disease: a clinical series. Cerebrovasc. Dis. 22, 109–115. doi: 10.1159/000093238

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D. Y., Son, J. P., Yeon, J. Y., Kim, G.-M., Kim, J.-S., Hong, S.-C., et al. (2017). Infarct pattern and collateral status in adult moyamoya disease. Stroke 48, 111–116. doi: 10.1161/STROKEAHA.116.014529

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, W. J., Jeong, S. K., Han, K. S., Lee, S. H., Ryu, Y. J., Sohn, C. H., et al. (2020). Impact of endothelial shear stress on the bilateral progression of unilateral moyamoya disease. Stroke 51, 775–783. doi: 10.1161/STROKEAHA.119.028117

PubMed Abstract | CrossRef Full Text | Google Scholar

Milewicz, D. M., Ostergaard, J. R., Ala-Kokko, L. M., Khan, N., Grange, D. K., Mendoza-Londono, R., et al. (2010). De novo ACTA2 mutation causes a novel syndrome of multisystemic smooth muscle dysfunction. Am. J. Med. Genet. A 152A, 2437–2443. doi: 10.1002/ajmg.a.33657

PubMed Abstract | CrossRef Full Text | Google Scholar

Starosolski, Z., Villamizar, C. A., Rendon, D., Paldino, M. J., Milewicz, D. M., Ghaghada, K. B., et al. (2015). Ultra high-resolution in vivo computed tomography imaging of mouse cerebrovasculature using a long circulating blood pool contrast agent. Sci. Rep. 5:10178. doi: 10.1038/srep10178

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsai, M., Kita, A., Leach, J., Rounsevell, R., Huang, J. N., Moake, J., et al. (2012). In vitro modeling of the microvascular occlusion and thrombosis that occur in hematologic diseases using microfluidic technology. J. Clin. Invest. 122, 408–418. doi: 10.1172/JCI58753

PubMed Abstract | CrossRef Full Text | Google Scholar

Urick, B. Y., Sanders, T. J., Hossain, S., Zhang, Y., and Hughes, T. J. (2019). Review of patient-specific vascular modeling: template-based isogeometric framework and the case for CAD. Arch. Comput. Methods Eng. 26, 381–404. doi: 10.1007/s11831-017-9246-z

CrossRef Full Text | Google Scholar

Wahlin, A., Ambarki, K., Birgander, R., Wieben, O., Johnson, K. M., Malm, J., et al. (2013). Measuring pulsatile flow in cerebral arteries using 4D phase-contrast MR imaging. AJNR Am. J. Neuroradiol. 34, 1740–1745. doi: 10.3174/ajnr.A3442

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeon, J. Y., Shin, H. J., Kong, D. S., Seol, H. J., Kim, J. S., Hong, S. C., et al. (2011). The prediction of contralateral progression in children and adolescents with unilateral moyamoya disease. Stroke 42, 2973–2976. doi: 10.1161/STROKEAHA.111.622522

PubMed Abstract | CrossRef Full Text | Google Scholar

Zipfel, G. J., Sagar, J., Miller, J. P., Videen, T. O., Grubb, R. L. Jr., Dacey, R. G. Jr., et al. (2009). Cerebral hemodynamics as a predictor of stroke in adult patients with moyamoya disease: a prospective observational study. Neurosurg. Focus. 26:E6. doi: 10.3171/2009.01.FOCUS08305

PubMed Abstract | CrossRef Full Text | Google Scholar