
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13098-8
10.1016/j.heliyon.2024.e37067
e37067
Research Article
Numerical analysis of the effect of Syringomyelia on cerebrospinal fluid dynamics
Edrisnia Hadis
Sharbatdar Mahkame m.sharbatdar@kntu.ac.ir
⁎
Faculty of Mechanical Engineering, K. N. Toosi University of Technology, Tehran, 19991-43344, Iran
⁎ Corresponding author. m.sharbatdar@kntu.ac.ir
29 8 2024
15 9 2024
29 8 2024
10 17 e3706722 6 2024
19 8 2024
27 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Spinal cord enlargement (SCE) includes conditions such as Syringomyelia, tumors, and tumor-like cases of demyelination, edema, or inflammation. These conditions involve fluid-filled cysts, known as syrinx, or masses of tissue, referred to as tumors, which cause increased pressure within the spinal cord (SC) and obstruct cerebrospinal fluid (CSF) circulation. To assess the impact of SCE location and diameter, we constructed fifteen computational SC models, each featuring a SCE placed in one of five probable locations with 20 %, 40 %, and 60 % stenosis. Our objective was to investigate how the location, diameter, and length of the SCE influence CSF velocity pattern and to identify the most critical location in the SC associated with this condition. The results indicated a velocity increase of 0.5 cm/(s) near models with 60 % stenosis. Importantly, SCE located from T1 to T5 exhibit a more pronounced reduction, exceeding 6.5, in the Womersley number. Our finding suggests that this region is the most vulnerable for SCE formation due to its significant impact on fluid circulation. The identification of specific locations within the SC associated with heightened risk can contribute to an improved understanding, treatment and management of SCE.

Keywords

Spinal cord enlargement (SCE)
Cerebrospinal fluid (CSF)
Computational fluid dynamic (CFD)
Spinal subarachnoid space (SSS)
Spinal cord (SC)
Tumor
==== Body
pmc1 Introduction

Spinal Cord Enlargement, SCE, diseases such as Syringomyelia, tumors, and tumor-like cases of demyelination, edema, or inflammation is a condition affecting the Spinal Cord, SC, characterized by the development fluid-filled cysts, known as syrinx, or tissue known as tumors within the cord tissue, respectively. Over time, these SCE expand, potentially causing damage to the SC tissue and leading to various neurological disorders. Symptoms may include weakness and decreased sensation. This condition is frequently associated with diseases that impede the flow of cerebrospinal fluid, CSF, within the Spinal Subarachnoid Space, SSS, such as SC injuries and Chiari type I malformation, CM1 [1]. Concerning physiological characteristics, the SSS surrounds the cord enclosed by the pia mater, a flexible membrane that permits the passage of CSF and tiny solutes. On the outer side, a rigid dura mater separates the subdural from the epidural space [2].

Recent studies have significantly advanced our understanding of CSF dynamics through various methodologies, generally falling into several categories. Firstly, research has explored the influence of anatomical features and physiological factors impacting CSF dynamics. Studies have shown how SC motion and the detailed structure of the SSS affect velocity fields and solute dispersion [3,4], and how micro-anatomical features in the cervical SSS improve local CSF velocity [5,6]. For example, Sass et al. [7] developed an anatomically accurate 3D model of the SSS with partially idealized nerve roots, and Khani et al. [8] demonstrated that nerve roots significantly alter the CSF's velocity field, steady streaming, and vortical features. Additionally, Yildiz et al. [9] showed that respiration and heart pulsations have varying effects on CSF velocity depending on the breathing patterns.

Secondly, CFD analyses have been utilized to model irregular flow patterns and evaluate the impact of anatomical variations [10,11]. However, comparative studies by Pahlavian et al. [12] and Yiallourou et al. [13] revealed that 4D PC-MRI measurements consistently have higher peak CSF velocities than CFD models. Montes et al. [14] introduced a time-averaged transport equation to address challenges in fluid-structure interaction (FSI) modeling, while Cardillo and Camporeale [15] constructed a comprehensive analytical FSI model revealing that deviations in Young's modulus values could increase CSF fluxes through the SC. Khani et al. [16] further validated CFD simulations and in vitro models in the context of intrathecal drug delivery.

Thirdly, the effect of spinal cord diseases, such as Chiari malformation type 1 (CM1), on CSF flow is studied. Research has highlighted elevated peak systolic and diastolic velocities in CM1 patients compared to healthy controls [17]. Furthermore, findings from Martin et al. [18] indicate robust inter-operator reliability for MRI-based CFD models of CSF flow in healthy subjects but diminished reliability in CM1 patients. Additionally, a new CSF hydrodynamic index based on maximum CSF pressure, derived from PC-MRI and CFD simulations, has proven effective in assessing CM1 patient conditions [19,20]. Tangen et al. [21] also emphasized the importance of body posture and drainage rates for effective blood clearance during lumbar drainage.

Table 1 summarizes the primary numerical studies on CSF dynamics related to Syringomyelia. Notably, the existing literature lacks detailed structures in the model. To the best of our knowledge, there has not been a comprehensive investigation into the influence of spinal cord enlargement size and location on CSF dynamics utilizing subject-specific spinal geometry. The principal aim of this study is to assess how these variables affect CSF hemodynamic parameters, utilizing computational fluid dynamics.Table 1 Summary of literature on CSF dynamics, including the region of interest, computational method used, geometry types, and key findings from each study.

Table 1Ref.	Full spine	FSI/CFD	SS/PSa	Method	Summary of results	
Bertram [22]	✓	One way
Finite element
ADINA R&D Inc.	Idealized (2D)	Thoracic syrinx measuring 140 mm in length - The degree of stenosis varied in specific areas, ranging from 0 % to 90 %.	The presence of bulk fluid flow adjacent to the syrinx led to substantial longitudinal pressure dissociation, often referred to as “suck".	
Cheng et al. [23]	–	CFD
Finite volume
ANSYS CFX	SS (3D)	Influence of arachnoiditis and its impact on CSF flow and pressure dynamics.	Arachnoiditis elevates flow resistance within the SSS, leading to a slight increase and shift in the CSF pressure drop across the affected region.	
Cheng et al. [24]	–	Two-way
ANSYS CFX and Mechanical	PS (3D)	Impact of FSI on CSF pressure and SC motion in both a normal and a stenosis model of the SSS.	FSI has a minimal impact on CSF pressure, with an approximately 1 % effect observed in both normal and stenosis models.	
Heil et al. [25]	✓	Two-way
Finite elements
MUMPS	Idealized (2D)	A poroelastic FSI model, encompassing the partially obstructed SSS, the syrinx, and the surrounding poroelastic tissues to study mass transfer dynamics between the SSS and the syrinx.	Poroelastic fluid-structure interaction can indeed lead to a slight increase in syrinx volume.	
Elliot et al. [26]	–	Two-way	idealized (2D)	Influence of Poisson's ratio on wave mode characteristics, analyzing dispersion behavior and mode forms across a wide wavelength spectrum.	Anticipated that with the ongoing radial expansion of the syrinx, the fluid sloshing within it and the resulting tension will diminish.	
Vinje et al. [27]	–	Two-way
Finite elements
FEniCS	Idealized (2D)	Evaluating Syrinx fluid velocity in a patient before and after craniovertebral decompression manipulating factors such as heart rates, syrinx diameters, and CSF flow velocities.	Syrinx fluid velocities are impacted by heart rate, CSF velocities, and syrinx diameter, but not by the frequency of oscillation.	
a Subject-Specific/Patient-Specific.

2 Method

Geometry and mesh. For this study, Sass et al.'s [7] publicly available, open-source 3D spine geometry has served as the foundation. For this purpose, A 23-year-old female was selected as a representative healthy subject to establish the CSF space geometry and flow parameters for the numerical model. She had no prior history of spinal deformities or CSF-related disorders. To obtain this geometry, high-resolution MR images are acquired to determine the entire intrathecal CSF space geometry, with a slice thickness of 1 mm, slice spacing of 0.499 mm, and in-plane voxel spacing of 0.547 × 0.547 mm2. The geometry is categorized into three sections as shown in Fig. 1, the first spans from FM to C2-C3, the second from C2-C3 to C7-T1, and the third from C7-T1 to T10-T11. The fluid domain is meshed using ANSYS meshing 20.1 (ANSYS Inc., Canonsburg, PA, USA). Three prismatic layers with the first layer thickness of 0.25 mm are used to increase the accuracy near the wall. then peak diastolic velocities are calculated and the current mesh model described in this context has demonstrated encouraging outcomes in regards to error percentage and computationally efficient processing time.Fig. 1 Schematic overview of the study: (a) BCs for three sections of the computational models, (b) 3D geometry of the fluid domain, (c) Detailed view of the fluid domain in the second section with SCE incorporation, and (d) Equation used to calculate the SCE diameter based on the stenosis percentage, transverse diameter of the dura and cord.

Fig. 1

Boundary conditions. Wall motion in the spinal canal induces fluid circulation; however, simulating fluid-structure interaction (FSI) and dynamic mesh requires significant computational resources. Given that FSI has a negligible impact on CSF dynamics [21], the dural and pial walls are modeled as rigid, with no-slip boundary conditions (BCs) applied. This simplification accurately represents the interaction between the CSF and canal walls. The base flow for each simulation is determined using MRI measurements and adjusted to ensure consistency across the geometry, including normal heart and respiration rates.

Domain division and flow specifications. Based on MRI measurements taken at three different locations for this subject-specific geometry, the fluid domain is divided into three distinct sections each simulated separately with mass flow rate, MFR, at the inlet and constant pressure at the outlet. The MFR for all sections is obtained from MRI measurements specific to this subject-specific geometry. Table 2 provides a comprehensive classification of the boundary conditions applied to each section.Table 2 Classification of sections for defining BCs in the computational models.

Table 2Section	region	Inlet	Outlet	Inlet BC	Outlet BC	
First	FM to C2-C3	C2-C3	FM	C2-C3 MFR	0 [Pa]	
Second	C2-C3 to C7-T1	C2-C3	C7-T1	C2-C3 MFR plus 0.0002679 [kg/s]	0 [Pa]	
Third	C7-T1 to T10-T11	C7-T1	T10-T11	C7-T1 MFR plus 0.00033 [kg/s]	1.06658 [Pa]	

For first section (FM to C2-C3) a MFR derived from MRI measurements [8] is applied at the inlet (C2-C3), with an outlet pressure of 0 Pa at the FM as shown in Fig. 1(a).

The second section (C2-C3 to C7-T1) incorporates an additional small MFR to account for wall movement. The combined flow rate from the first section is applied at the C2-C3 inlet, with the pressure at the C7-T1 outlet set to 0 Pa.

Finally, for the third section, an MRI-derived MFR is applied at the C7-T1 inlet, with a constant positive pressure at the T10-T11 outlet to align the peak CSF velocities observed along the geometry [8].

Computational model. The continuity and Navier-Stokes equations (Eqs. (1), (2), respectively) are solved by CFD software ANSYS-CFX 20.1 (ANSYS Inc., Canonsburg, PA, USA) as provided below:(1) ∇.u→=0

(2) ρ∂u→∂t+ρ(u→.∇)u→+∇p−μΔu→=0

where u→ is the fluid velocity vector and μ,ρ and p represent the local dynamic viscosity, density, and pressure, respectively. The timestep is assumed 0.005 s and maximum residual mean square errors are set to 1 × 10−4. Additionally, CSF is considered an incompressible Newtonian fluid with a viscosity of 0.6913 mPa s and a density of 993.3 kg/m3 at body temperature, the same as water [8,28]. To ensure periodic repetitive states, two cardiac cycles are conducted, and the peak diastolic velocity is compared between them. The negligible 2.88 % error in Figure S 1 shows that the use of the first cycle is justifiable.

SCE modeling. In total, 15 models are created, distributed across three distinct sections, three in the first section from C1 to C3, six in the second section, C3-C4 to C6 and C4-C5 to C7, and six in the third section, T1 to T3 and T3 to T5, shown in Fig. 2. The almond shaped SCE as a separated STL file, exhibiting stenosis percentages of 20, 40, and 60 percent, are meticulously inserted and merged with cord at various points along the SC, incorporating differing lengths to ensure a comprehensive exploration. This computational approach aims to encompass a spectrum of scenarios reflecting the diverse manifestations of SCE disease, with the selection of positions roughly aligned with those detailed in Ref. [29] based on morphological evaluation of the SCE. In Fig. 1(c), a representative model of a SCE integrated into the second section is visually depicted, offering a tangible illustration of this incorporation process. Moreover, the formulation for determining the maximum transverse diameter of SCEs, derived from the interplay of stenosis percentage, the maximum transverse diameters of dura, and cord cross-sections, is thoughtfully elucidated in Fig. 1(d). This formula is adapted from Edrisnia et al. [30] who presented it for calculating stenosis percentage in coronary arteries.Fig. 2 Comparative Analysis of SCE Incorporation across stenosis regions specified by dash line squares, showcasing three degrees of stenosis in each column. (a) general schematic of the three sections, (b) five locations with 20 % of stenosis, (c) five locations with 40 % of stenosis, (d) five locations with 60 % of stenosis.

Fig. 2

hydrodynamic parameters. The spine's geometric and hydrodynamic parameters are assessed at 5 mm intervals. The SSS cross-sectional area, Acs, is computed as the difference between the cross-sectional area of the dura, Ad, and the cord, Ac. The hydraulic diameter, DH=4AcsPcs, for the internal flow within the tube, is calculated based on the cross-sectional area and wetted perimeter, Pcs=Pd+Pc. All these parameters are determined using a user-defined function compiled in ANSYS CFD-Post (ANSYS Inc., Canonsburg, PA, USA). The Reynolds number, Re, is based on the peak flow rate, while the Womersley number, α, relies on the hydraulic diameter. Reynolds number is calculated as Re=QsysDHνAcs, where Qsys represents the temporal maximum of the local flow at each axial interval along the spine, and ν is the fluid's kinematic viscosity. To assess the presence of laminar flow, the Reynolds number is evaluated based on previous research (Re<2300) for peak systolic flow along the spine [7]. The Womersley number is measured to determine the ratio of unsteady inertial forces to viscous forces. It is calculated as α=Dh2ω/ν, where ω is the angular velocity of the volume flow waveform, ω=2π/T. A Womersley number greater than 5 for oscillatory flows indicates a transition from parabolic to “m-shaped” velocity profiles. In our study, Womersley number is in the range of 5–20.

3 Result

Grid study and mesh generation. A grid analysis is performed for the fluid domain within the stenosis model, featuring a 60 % stenosis in the second section. Mesh generation is carried out using ANSYS-Meshing 20.1 (ANSYS Inc., PA, USA). The mesh configuration involved a uniform tetrahedral mesh for the CSF flow analysis. Subsequently, identical element sizes and layer quantities are applied to the remaining stenotic geometries. The mesh was refined by applying smaller elements near the dura and spinal cord walls due to capture the significant flow gradients in these regions. To ensure mesh independence, a grid sensitivity analysis is conducted by comparing simulations with varying cell numbers. The analysis involved calculating the average velocity at the cross-section of the stenosis midpoint for different mesh types. Table 3 illustrates that the medium mesh, with a node number of approximately 36,000 and cell number of 125,000, provided results with less than 0.2 % error relative to the finest mesh, confirming its adequacy for accurate simulations.Table 3 Numerical sensitivity study, showing the mesh type data and difference between average velocity calculated at the midpoint of the stenosis for the second section.

Table 3Mesh type	Cell number	Node number	First layer thickness (mm)	Prismatic layer	error (%)	
Coarse	95,000	25,000	0.3	2	1.51	–	
Medium	125,000	36,000	0.25	3	0.2	
Fine	490,000	142,000	0.1	4	–	

Verification. For verification purposes, the peak diastole velocity in the z-direction is compared between the computed values obtained from the current simulation and the data from earlier investigations [8], as illustrated in Fig. 3. The comparison involved assessing the peak diastole CSF velocity in z-direction throughout the whole geometry in 1 [mm] intervals. The total average error is less than %14.Fig. 3 Verification of the results. Computed peak diastolic CSF velocity in the z direction throughout the entire geometry compared with data of Khani et al. [8].

Fig. 3

Velocity distribution. In Fig. 4, the distributions of velocity magnitude during peak diastole are depicted for nine distinct models, each featuring SCE-induced stenosis across three specific sections: C1 to C3, C4-C5 to C7, and T3 to T5. The severity of stenosis varies among these models, categorized as 20 %, 40 %, and 60 %. There is a nearly 50 % decline in the lower path in the 60 % of stenosis in Fig. 4(c) compared to the 20 % stenosis in Fig. 4(a) which means the pressure increases in this region An increase in a jet-shaped velocity profile is observed in the upper path of Fig. 4(c). Furthermore, as the stenosis severity increases from %20 (Fig. 4(d)) to %40 (Fig. 4(e)), there is a lateral decrease in velocity magnitude, resulting to a substantial reduction of 70 %, ultimately reaching zero velocity. However, at both the inlet and outlet, all cases display consistent trends and values due to consistent BCs. Moreover, in Fig. 4(i), where there is 60 % stenosis in T3 to T5 section, the velocity in the lower path decreases by 25 % compared to Fig. 4(g).Fig. 4 Velocity magnitude distributions at peak diastole in nine distinct models, each containing three SCE regions across three different sections, C1 to C3 (a, b, and c), C4-C5 to C7 (d, e, and f), and T3 to T5 (g, h, and i). Stenosis percentages vary, with (c, f, and i) representing 20 % stenosis, (b, e, and h) 40 % stenosis, (a, d, and g) 60 % stenosis.

Fig. 4

Comparative Analysis of Systolic and Diastolic Velocities.Table 4 provides information about the peak diastolic and systolic velocity magnitude both exactly above and below each SCE in the locations specified in the table. Generally, an increase in velocity is observed near the stenotic regions as the SCE expands from 20 % to 60 % of stenosis. These changes are particularly significant in the third section in peak diastole compared to other areas. For instance, in the first section, C1 to C3, the SCE causes an approximately 0.1 cm/(s) increase in systolic and diastolic velocity magnitude above the SCE when the stenosis area triples. In the second region, C3-C4 to C6 and C4-C5 to C7, velocity magnitudes show less than a 0.1 cm/(s) change. The model without any SCE shows higher velocity in both systole and diastole peaks in the same location as below the SCE in C4-C5 to C7. However, in the third area with narrowing in T1 to T3 and T3 to T5, an expansion of the SCE from 20 % to 60 % results in a growth of about 0.2–0.3 cm/(s). Furthermore, in peak diastole the existence of SCE increases the velocity by almost 0.5 cm/(s) both above and below the SCE, while, in peak systole the existence of SCE decreases the velocity by about 0.05–0.3 cm/(s).Table 4 Peak velocity magnitude above and below each SCE, with A and B representing above and below the SCE, respectively. All velocities are expressed in [cm/s].

Table 4Location
Stenosis	C1 to C3	C3-C4 to C6	C4-C5 to C7	T1 to T3	T3 to T5	
A	B	A	B	A	B	A	B	A	B	
Peak diastolic velocity	
0 %	0.597	0.817	0.935	1.146	1.103	1.485	1.857	1.523	1.476	1.258	
20 %	0.596	0.804	0.939	1.153	1.099	1.419	1.821	1.525	1.464	1.248	
40 %	0.625	0.819	0.952	1.182	1.089	1.382	1.860	1.615	1.523	1.308	
60 %	0.692	0.867	0.976	1.239	1.108	1.426	1.918	1.854	1.728	1.443	
Peak systolic velocity	
0 %	1.040	1.480	1.293	1.598	1.516	2.009	1.853	1.503	1.482	1.236	
20 %	1.022	1.424	1.299	1.608	1.511	1.902	1.790	1.503	1.460	1.226	
40 %	1.061	1.401	1.312	1.645	1.502	1.866	1.810	1.576	1.494	1.297	
60 %	1.150	1.486	1.340	1.718	1.521	1.951	1.856	1.767	1.688	1.401	

Fig. 5 illustrates the velocity distribution at two time steps for two distinct geometries: one without stenosis (Fig. 5(a and b)) and another with 60 % stenosis (Fig. 5(c and d)). Generally, it can be seen that in the presence of 60 % of stenosis during both systole and diastole (Fig. 5(c and d)), there is no significant impact on the low-velocity pattern in the lower path. However, in the upper path, the velocity increases by 25 % during systole compared to diastole. Moreover, at both time steps, the presence of 60 % of stenosis results in a reduction of velocity by more than 50 % in the lower path in comparison to the without stenosis model.Fig. 5 Comparative analysis of systolic and diastolic velocities for a single location in C3-C4 to C6. (a) without stenosis in systole, (b) without stenosis in diastole, (c) with 60 % stenosis in systole, and (d) with 60 % stenosis in diastole.

Fig. 5

Hydrodynamic parameters.Fig. 6 presents hydrodynamic results for a total of seven cases. Fig. 6(b) and (c) illustrate three distinct models featuring unique regions for SCE with varying severities. In addition, five different potential locations [29] all with the same 60 % degree of stenosis are depicted in Fig. 6(d) and (e). In a known position spanning from C4-C5 to C7 shown in Fig. 6(b), the model with 60 % stenosis exhibits a notable lower local minimum Womersley number as 7.68 compared to the corresponding value for 20 % stenosis accounted for 10.78. In addition, for all cases featuring SCEs at distinct locations with the same 60 % stenosis, the Womersley number ranges from α = 23.16 to α = 4.5, as depicted in Fig. 6(d), with a local minimum occurring in the SCE region. Notably, the presence of the SCE from T1 to T3 has resulted in a more pronounced reduction in the Womersley number compared to the other models shown in Fig. 6(d), with a reduction of more than 6.5 units in the same location. Concerning the Reynolds number, Fig. 6(b) illustrates the Reynolds number in peak systole (t = 0.6 s) for three various degrees of stenosis at a specific location. The noticeable point observed in the line graph is that the three models do not exhibit identical values, except for those near the inlet and outlet, owing to their similar BCs as a consequence of the Reynolds number above and below the stenosis. Furthermore, an increase in stenosis percentage induces greater variability in Reynolds values, resulting in nearly a tenfold increase in the minimum (Re = 120) and maximum (Re = 150) values compared to the 20 % stenosis scenario. These changes in velocity patterns and values become even more pronounced as the SCE diameter increases. Regarding Fig. 6(e), depicting 60 % stenosis in five distinct locations, the most significant change was observed from T1 to T3. This outcome was anticipated since the maximum alteration in the Womersley number occurred for this particular model in comparison to the other areas. To provide precise details, the graph reached its highest point at Re = 217, situated between T1 and T3 where the SCE was positioned, 15 cm below the region of FM. Moreover, Fig. 6(f, g, and h) compares the systolic and diastolic Reynolds number at three distinct locations across all three sections. The comparative analysis shows that the stenosis in C1 to C3 and C4-C5 to C7 has a more pronounced effect, almost twice as much, on the Reynolds number trend during systole relative to diastole in Fig. 6(f and g). However, this effect is less significant in the third section, as illustrated in Fig. 6(h).Fig. 6 Distribution of hydrodynamic parameters along the spine: (a) 3D geometry with axial plane locations of added SCEs, (b, c) Womersley and Reynolds numbers for Three stenosis percentages within C4-C5 to C7 in peak diastole, (d, e) Womersley and Reynolds numbers for 60 % stenosis at five different spine positions in peak diastole, (f) Systolic and diastolic Reynolds number comparison at C1 to C3, (g) Systolic and diastolic Reynolds number comparison at C4-C5 to C7, (h) Systolic and diastolic Reynolds number comparison at T3 to T5.

Fig. 6

Pressure distribution.Fig. 7 demonstrates the pressure distribution for two models, Fig. 7(a) shows section three without stenosis and Fig. 7(b) shows section three with 60 % of stenosis in T1 to T3. What stands out from this comparative analysis is that the presence of the stenosis reduces the pressure value both before and after the stenosis.Fig. 7 Pressure distribution at peak systole: (a) Section three without stenosis and (b) Section three with 60 % of stenosis in T1 to T3.

Fig. 7

4 Discussion

The general goal of this work is to investigate the effect of SCE on CSF motion in SSS by adding SCEs to the subject-specific geometry. This study can be considered as a reproduction of the work presented by Khani et al. [8], with modifications to the geometry to include SCEs and substantial simplifications in the modelling approach, such as applying BCs. SCEs are added to potential locations from the cervical to the thoracic segment of the SC [1,29]. The results have been studied over a range of model geometries with the CSF given correct properties [8]. This study does not consider the poroelastic property of the cord, therefore the changes in SCE volume due to FSI and poroelasticity of the cord studied in Refs. [1,22,25] is not considered. However, it has shown that the change in SCE volume in the full FSI model is very small [25]. In Fig. 7, the presence of the SCE increases the pressure drop, and the magnitude of pressure after the SCE in Fig. 7(b) is almost halved compared to Fig. 7(a). This finding is supported by previous studies that test the presence of microstructures [22,23]. Moreover, the peak CSF velocity in Fig. 5 depicts that flow patterns are more complex in geometry with SCE and velocity is higher adjacent to the SCE as showed by existing MRI studies [31]. As the effects of FSI on CSF pressure are shown negligible in stenosis models in previous studies [24], the walls are considered rigid. Despite a substantial SCE, the outer cord diameter was maintained unchanged locally in previous studies [1,27] and they did not consider the effect of SCE on CSF motion in the SSS. For the first time in the literature, our simulations show that the velocity in systole and diastole and Reynolds number in SSS are affected with the addition of SCE. In general, the measured values for velocity and Reynolds number in Table 4 and Fig. 6 were comparable with those measured in vivo in humans [8].

Limitation and Future Remarks. The primary limitation of our study is that, for clinical application, the geometry reconstruction should be done automatically using deep learning algorithms rather than manual segmentation as in previous studies [7,8]. Secondly, the size of the SCEs added to the spinal cord is limited to 20 %, 40 %, and 60 %. Therefore, future studies should consider a wider range of sizes. Thirdly, microstructures such as the denticulate ligaments, trabeculae, and nerve roots were excluded because including them in all models was not computationally efficient. However, it has been shown that the presence of these microstructures does not significantly change the velocity profiles. Fourthly, our models did not account for cord motion and permeability. Instead, fluid circulation was simulated by MFR and pressure BCs. Moreover, Recent work by Versaci et al. [32] illustrates how fuzzy similarity techniques can enhance anomaly detection and data analysis in complex datasets, offering valuable insights for improving CFD methodologies. Therefore, the effect of this technique in CSF flow analysis can be studied in future. Finally, a comparison between our approach and computational-clinical PC-MRI studies on cases with the same disease should be conducted for direct validation.

5 Conclusion

In this study, we proposed a proof-of-concept framework to assess the impact of SCE, which is a neurological disorder characterized by fluid-filled cysts, on CSF dynamics. To achieve this goal, we developed fifteen computational models, each featuring SCE placement in one of five locations and with mild, moderate, and severe stenosis. The aim was to investigate how SCE presence influences CSF velocity pattern and find the most threatening location in the SC for this disease. In summary, the identified dangerous locations that have been located in T1 to T5 within the SC provide a foundation for targeted interventions and further research in this field. However, this research has some limitations that can be addressed in future studies such as considering cord motion and permeability and including fine structures in the canal to find the pattern of fluid circulation more accurately. Moreover, the geometry of disease may be studied to have the real shape of SCE instead of adding the idealized model to a healthy SC.

Data availability

Data will be made available on request.

CRediT authorship contribution statement

Hadis Edrisnia: Writing – original draft, Validation, Methodology, Formal analysis, Data curation. Mahkame Sharbatdar: Writing – review & editing, Supervision, Methodology, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Supplementary data

Figure S1peak CSF velocity comparison between first and second cycle.

Appendix A Supplementary data

The following is the Supplementary data to this article:Multimedia component 1

Multimedia component 1

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e37067.
==== Refs
References

1 Bertram C.D. Heil M. A poroelastic fluid/structure-interaction model of cerebrospinal fluid dynamics in the cord with syringomyelia and adjacent subarachnoid-space stenosis J. Biomech. Eng. 139 2017 10.1115/1.4034657
2 Reina M.A. Collier C.B. Prats-Galino A. Puigdellívol-Sánchez A. MacHés F. De Andrés J.A. Unintentional subdural placement of epidural catheters during attempted epidural anesthesia: an anatomic study of spinal subdural compartment Reg. Anesth. Pain Med. 36 2011 537 541 10.1097/AAP.0b013e31822e0e8c 21989150
3 Loth F. Yardimci M.A. Alperin N. Hydrodynamic modeling of cerebrospinal fluid motion within the spinal cavity J. Biomech. Eng. 123 2001 71 79 10.1115/1.1336144 11277305
4 Stockman H.W. Effect of anatomical fine structure on the dispersion of solutes in the spinal subarachnoid space J. Biomech. Eng. 129 2007 666 675 10.1115/1.2768112 17887892
5 Pahlavian S.H. Yiallourou T. Tubbs R.S. Bunck A.C. Loth F. Goodin M. Raisee M. Martin B.A. The impact of spinal cord nerve roots and denticulate ligaments on cerebrospinal fluid dynamics in the cervical spine PLoS One 9 2014 1 14 10.1371/journal.pone.0091888
6 Tangen K.M. Hsu Y. Zhu D.C. Linninger A.A. CNS wide simulation of flow resistance and drug transport due to spinal microanatomy J. Biomech. 48 2015 2144 2154 10.1016/j.jbiomech.2015.02.018 25888012
7 Sass L.R. Khani M. Natividad G.C. Tubbs R.S. Baledent O. Martin B.A. A 3D subject-specific model of the spinal subarachnoid space with anatomically realistic ventral and dorsal spinal cord nerve rootlets Fluids Barriers CNS 14 2017 1 16 10.1186/s12987-017-0085-y 28132644
8 Khani M. Sass L.R. Xing T. Keith Sharp M. Balédent O. Martin B.A. Anthropomorphic model of intrathecal cerebrospinal fluid dynamics within the spinal subarachnoid space: spinal cord nerve roots increase steady-streaming J. Biomech. Eng. 140 2018 10.1115/1.4040401
9 Yildiz S. Thyagaraj S. Jin N. Zhong X. Heidari Pahlavian S. Martin B.A. Loth F. Oshinski J. Sabra K.G. Quantifying the influence of respiration and cardiac pulsations on cerebrospinal fluid dynamics using real-time phase-contrast MRI J. Magn. Reson. Imag. 46 2017 431 439 10.1002/jmri.25591
10 Roldan A. Wieben O. Haughton V. Osswald T. Chesler N. Characterization of CSF hydrodynamics in the presence and absence of tonsillar ectopia by means of computational flow analysis Am. J. Neuroradiol. 30 2009 941 946 10.3174/ajnr.A1489 19299486
11 Linge S.O. Haughton V. Løvgren A.E. Mardal K.A. Langtangen H.P. CSF flow dynamics at the craniovertebral junction studied with an idealized model of the subarachnoid space and computational flow analysis Am. J. Neuroradiol. 31 2010 185 192 10.3174/ajnr.A1766 19729542
12 Pahlavian S.H. Bunck A.C. Loth F. Tubbs R.S. Yiallourou T. Kroeger J.R. Heindel W. Martin B.A. Characterization of the discrepancies between four-dimensional phase-contrast magnetic resonance imaging and in-silico simulations of cerebrospinal fluid dynamics J. Biomech. Eng. 137 2015 10.1115/1.4029699
13 Yiallourou T.I. Kröger J.R. Stergiopulos N. Maintz D. Martin B.A. Bunck A.C. Comparison of 4D phase-contrast MRI flow measurements to computational fluid dynamics simulations of cerebrospinal fluid motion in the cervical spine PLoS One 7 2012 10.1371/journal.pone.0052284
14 Gutiérrez-Montes C. Coenen W. Lawrence J.J. Martínez-Bazán C. Sánchez A.L. Lasheras J.C. Modelling and direct numerical simulation of flow and solute dispersion in the spinal subarachnoid space Appl. Math. Model. 94 2021 516 533 10.1016/j.apm.2021.01.037
15 Cardillo G. Camporeale C. Modeling fluid–structure interactions between cerebro-spinal fluid and the spinal cord J. Fluid Struct. 102 2021 103251 10.1016/j.jfluidstructs.2021.103251
16 Khani M. Burla G.K.R. Sass L.R. Arters O.N. Xing T. Wu H. Martin B.A. Human in silico trials for parametric computational fluid dynamics investigation of cerebrospinal fluid drug delivery: impact of injection location, injection protocol Physiology 2022 10.1186/s12987-022-00304-4
17 Rutkowska G. Haughton V. Linge S. Mardal K.A. Patient-specific 3D simulation of cyclic CSF flow at the craniocervical region Am. J. Neuroradiol. 33 2012 1756 1762 10.3174/ajnr.A3047 22517282
18 Martin B.A. Yiallourou T.I. Pahlavian S.H. Thyagaraj S. Bunck A.C. Loth F. Sheffer D.B. Kröger J.R. Stergiopulos N. Inter-operator reliability of magnetic resonance image-based computational fluid dynamics prediction of cerebrospinal fluid motion in the cervical spine Ann. Biomed. Eng. 44 2016 1524 1537 10.1007/s10439-015-1449-6 26446009
19 Gholampour S. Gholampour H. Correlation of a new hydrodynamic index with other effective indexes in Chiari I malformation patients with different associations Sci. Rep. 10 2020 1 13 10.1038/s41598-020-72961-0 31913322
20 Pahlavian S.H. Loth F. Luciano M. Oshinski J. Martin B.A. Neural tissue motion impacts cerebrospinal fluid dynamics at the cervical medullary junction: a patient-specific moving-boundary computational model Ann. Biomed. Eng. 43 2015 2911 2923 10.1007/S10439-015-1355-Y 26108203
21 Tangen K. Narasimhan N.S. Sierzega K. Preden T. Alaraj A. Linninger A.A. Clearance of subarachnoid hemorrhage from the cerebrospinal fluid in computational and in vitro models Ann. Biomed. Eng. 44 2016 3478 3494 10.1007/s10439-016-1681-8 27384938
22 Bertram C.D. Evaluation by fluid/structure-interaction spinal-cord simulation of the effects of subarachnoid-space stenosis on an adjacent syrinx J. Biomech. Eng. 132 2010 1 15 10.1115/1.4001165
23 Cheng S. Stoodley M.A. Wong J. Hemley S. Fletcher D.F. Bilston L.E. The presence of arachnoiditis affects the characteristics of CSF flow in the spinal subarachnoid space: a modelling study J. Biomech. 45 2012 1186 1191 10.1016/j.jbiomech.2012.01.050 22386041
24 Cheng S. Fletcher D. Hemley S. Stoodley M. Bilston L. Effects of fluid structure interaction in a three dimensional model of the spinal subarachnoid space J. Biomech. 47 2014 2826 2830 10.1016/j.jbiomech.2014.04.027 25005435
25 Heil M. Bertram C.D. A poroelastic fluid-structure interaction model of syringomyelia J. Fluid Mech. 809 2016 360 389 10.1017/jfm.2016.669
26 Elliott N.S.J. Lucey A.D. Lockerby D.A. Brodbelt A.R. Fluid–structure interactions in a cylindrical layered wave guide with application in the spinal column to syringomyelia J. Fluid Struct. 70 2017 464 499 10.1016/j.jfluidstructs.2016.11.007
27 Vinje V. Brucker J. Rognes M.E. Mardal K.A. Haughton V. Fluid dynamics in syringomyelia cavities: effects of heart rate, CSF velocity, CSF velocity waveform and craniovertebral decompression NeuroRadiol. J. 31 2018 482 489 10.1177/1971400918795482 30114970
28 Gupta S. Soellinger M. Grzybowski D.M. Boesiger P. Biddiscombe J. Poulikakos D. Kurtcuoglu V. Cerebrospinal fluid dynamics in the human cranial subarachnoid space: an overlooked mediator of cerebral disease. I. Computational model J. R. Soc. Interface 7 2010 1195 1204 10.1098/rsif.2010.0033 20236960
29 Yeo X.J. Cheng X.S. Hemley X.S. Lee X.B.B. Stoodley X.M. Bilston X.L. Characteristics of CSF Velocity-Time Profile in 2017
30 Edrisnia H. Sarkhosh M.H. Mohebbi B. Parhizgar S.E. Alimohammadi M. Non-invasive fractional flow reserve estimation in coronary arteries using angiographic images Sci. Rep. 14 2024 1 15 10.1038/s41598-024-65626-9 38167627
31 Struck A. Radiology V.H.- Undefined 2009, Idiopathic syringomyelia: phase-contrast MR of cerebrospinal fluid flow dynamics at level of foramen magnum Pubs.Rsna.OrgAF Struck, VM HaughtonRadiology, 2009•pubs.Rsna.Org 253 2009 184 190 10.1148/radiol.2531082135
32 Versaci Mario A fuzzy similarity-based approach to classify numerically simulated and experimentally detected carbon fiber-reinforced polymer plate defects Sensors 22 11 2022 4232 35684853
