Computational Modeling and Simulation of Stenosis of Aqueduct of Sylvius Due to Brain Tumor

1 Background: Stenosis of cerebral aqueduct (CA) is featured in many studies related to elevated 2 intracranial cerebral pressures (ICP). It also presents a challenging situation to clinicians. Compressive 3 forces play a lead role in pathological situations like tumor presence and hence can cause obstruction to 4 the flow of cerebrospinal fluid (CSF). Due to this barrier, excessive retention of CSF in ventricles can 5 occur. This in turn could contribute to increased pressure gradients inside the cranium. In literature, most 6 of the numerical models are restricted to modeling the CSF flow by considering ventricle walls as rigid 7 material unlike its behavior a deformable character. This paper, therefore, addresses the same from a 8 holistic perspective by taking into consideration the dynamics of the flexible character of the ventricular 9 wall. This adds to the novelty of this work by reconstructing an anatomically realistic ventricular wall 10 behavior. To do this, the authors aim to develop a computational model of stenosis of CA due to brain 11 tumor by invoking a fluid-structure interaction (FSI) method. The proposed 3D FSI model is simulated 12 under two cases. First, simulation of pre-stenosis case with no interaction of tumor forces and secondly, a 13 stenosis condition together-with dynamic interaction of tumor forces. 14 Results: Comparing the forces with and without tumor reveals a marked obstruction of CSF outflow post 15 third ventricle and the cerebral aqueduct. Not only this but a drastic rise of CSF velocity from 21.2 mm/s 16 in pre-stenosis case to 54.1 mm/s stenosis case is also observed along with a net deformation increase of 17 0.144 mm on walls of ventricle. 18 Conclusions: This is a significant contribution to brain simulation studies for pressure calculations, 19 wherein the presence of tumors is a major concern.


BACKGROUND 1
Brain tumors can be of benign or malignant. Compressive forces of brain tumor may constrict the 2 flow of Cerebrospinal fluid (CSF), thereby causing stenosis of cerebral aqueduct (CA). Likewise, during 3 stenosis of CA, obstructive hydrocephalus (which is one of the direct consequences of constriction caused 4 by brain tumor on walls of CA) could be seen as well. Understanding the underlying mechanisms of 5 stenosis of CA, nature, pathophysiology and biomechanics of brain tumor and obstructive hydrocephalus 6 viz a viz their relationship with each other can be understood well. The study of stenosis of CA can also 7 aid in better understanding the rise of intracranial pressures (ICP) under special circumstances. According 8 to a survey conducted by [1] around 23,890 adults in the United States of America are diagnosed with 9 brain cancers this year; and whereas, every year about a million Americans are affected by Hydrocephalus 10 [2]. Stenosis of cerebral aqueduct (CA) is an important domain and catchy topic clinically to discuss. 11 Since continuous monitoring of CSF pressure (and ICP) inside the cranium usually involves surgical 12 interventions. The invasive mechanisms for clinically sensitive procedures are risky and cannot be used, 13 as a matter of routine. However, a numerical-based model may address this clinical and practical 14 limitation in a non-invasive manner. This can provide clinicians better analytics methods to understand 15 brain tumor interaction with walls of CA and its effects, including stenosis without surgical interventions. 16 CSF modeling exhibits fluid flow in brain noninvasively. For this purpose, computational fluid dynamics 17 (CFD) has been widely used to model fluid flows for biomedical applications. Many studies provide 18 evidence for the same including Jacobson et al. [3] which consists of a simple cylinder to study the flow 19 of CSF through of CA. This study suggests that pressure difference of at least 1.1 Pa is required to drive 20 the CSF from CA. In another study [4], a synchronized framework is proposed by considering two 21 models, a cylindrical rigid wall model and an elastic wall model segmented from MRI data to study flow 22 of CSF in CA. The spatial domain was digitized using Immersed Boundary method (IBM) [5]. Pulsatile  A. Linninger et al. [6] proposed a 2D rigid model of human ventricular system (HVS) where a pressure 3 drop of less than 2 Pa and velocity of 7.3 mm/s is reported. In one of the extensive studies on CSF flow 4 employing HVS, Howden et al. [7] 3D Ventricles were segmented from MRI data and CSF was treated as 5 incompressible Newtonian fluid [11]. Pulsatile velocity inlet and zero-gauge pressure at outlet boundary 6 conditions were used. The maximum velocity and Reynolds number was found to be 11.38 mm/s and 15 7 in CA respectively. Similarly, in a study by Vartan Kurtcuoglu et al. [8]., authors provided a guided 8 framework using simplified geometry of ventricles and compared their case with the simplified geometry 9 of stenosed. This study concluded a rift in pressure amplitudes and peak velocities change relative to the 10 case of stenosis. The study further concludes that under a stenosed CA, pressure amplitude in the lateral 11 ventricles increase as stenosis level increase. Study of Jacobson et al. [9] also suggests that under mild 12 stenosis, pressure required to initiate CSF flow increases to 125 Pa. A conclusion of this work reveals an 13 increase in the "Transmantle pressure" pivotally, even with little changes in intraventricular pressures. 14 One such approach is Fluid structure interaction (FSI

CSF Flow Field and Reynold's Number under Stenosis of CA 7
Velocity field is one of the important parameters in CSF modeling, CSF peak velocity always occurs 8 in the cross-section of CA due to narrow pathway. The simulation results are presented for two cases. 9 First case pertains to FSI simulation of ventricular body with no interaction of tumor forces. corresponds to the maximum stenosis of CA due to which velocity rises drastically to 54.1 mm/s 16 depicting that the duct has deformed considerably, relative to the previous stenosis levels. Fig 7 (d), 17 shows that the velocity vectors at the point constriction have started to change direction as well. This 18 suggests that the deformable property of ventricle should not be overlooked because it has visible effect 19 on the biomechanics of the system. 20 In addition to this, another useful parameter to gauge the dynamics of CSF under the scenario of stenosis 1 of CA was to find out Reynold's number of the flow. Reynold's number is another important parameter 2 which signifies the transition of flow from laminar to turbulent. It is given by: a healthy brain where no tumor is present. In this case, a maximum pressure of 2.5 Pa is observed. It is 8 also observed that the pressure variations in large cavities such as lateral ventricles were found to be 9 spatially uniform. This is perhaps because of their large size. Pressure drop of 0.2 Pa and 1.7 Pa is found 10 across the third ventricle and CA. A large pressure drop across CA is observed to due narrow pathway 11 resulting in greater pressure drop. 12

Deformation on Walls of Ventricles. 5
Since walls of ventricles have been modeled as elastic therefore their deformable character results in 6 deformation. Figure 9 show the deformation field obtained on the walls of ventricles under two cases. Case 2 : Fig 9 (b) shows a tumor specific case. In this particular case, a drastic increment in deformation 1 of 0.146 mm in lateral ventricles is observed. 2 This deformation shows that once the CSF flow is obstructed considerably, pressure in lateral ventricles 3 increase and deform the lateral ventricles. It also creates a requisite compliance for the excessive retention 4 of CSF fluid which cannot flow towards fourth ventricles due to constriction. Since this particular case 5 has not be modeled even in 2D hence the values presented herein can only be interpreted intuitively.

Effect of Considering CSF Pulsatility on Stenosis of CA 9
The pulsatile component of CSF is induced due to the blood flow in cerebral vasculature during a 10 cardiac cycle. Pulsatile mass flow inlet is provided using equation (10). Figure 10 shows the simulation 11 results because of cardiac cycle on the flow of CSF in ventricles. 12

DISCUSSION 5
The proposed study presents a 3D FSI model of stenosis of CA due to the compressive forces of 6 brain tumor. Our work is an effort to provide a basic framework from which any ventricle-CSF 7 interaction can be studied. For better visualization, two cases are simulated, with and without interaction 8 of tumor. This comparison allows us to visualize differences between the two conditions. For these two It is important to discuss here how the results were validated. Validation of results were carried out via 3 two independent sources. First, we validated our results with the past models of FSI presenting case 4 studies of CSF flow parameters. Second source of validation was PC-MRI studies, taken from the 5 literature, which provide experimental results for CSF flow parameters. Table 3 provides the results of 6 present case and compares it with previous FSI-based papers and Table 4  al. [7] are limited to-date as they consider ventricular surface as rigid. Table 4 presents some basic Cine- From a theoretical and practical perspective, we conclude that various flow parameters, magnitude and 21 location of tumor forces affect the outcome. One important parameter is Transmantle pressure. Higher the 22 Transmantle pressure higher would be the distention/deformation on the ventricular body. The same has 23 also been alluded in the study by Holmlund et al. [23][24]. Whereas, CSF pulsatility amplifies and 24 increases Transmantle pressures. Hence another point of conclusion is the physiological role and 1 importance of pulsatility in 3D FSI calculation of CSF in intraventricular and intracranial regime.

LIMITATIONS 2
There are limitations associated with this study too, which needs to be discussed here. First limitation 3 is the factor of gravity. Gravity affects the flow of fluids generally. However, if the flow is single-phase 4 and fluid is incompressible (which CSF is in proposed case) then gravity produces little effect on the flow 5 field because when there is no net change in density, the fluid forces will be exactly the same and counter-6 acting against the gravity forces. 7 Second is the factor of slight pulsations caused by respiration (similar to cardiac induced pulsations). In 8 cases where breathing is normal, CSF movement is little influenced by normal thoracic breathing [25], 9 however, if deep excursions and abdominal breathing is considered then it can be presumed that it may 10 induce some effect on the CSF flow [26]. 11 Third factor which was in the larger perspective but was not taken into account in this study is that we 12 have mainly focused on the tumor forces due to the solid core contributing body forces and have 13 neglected the fluid component. Though tumor 's solid core is essentially made up of extracellular matrix 14 and collagen fibers which by far has the largest contribution towards tumor forces and can be assumed to 15 approximate around this force, however, taking fluid forces also may introduce further improvement in 16 the model. It also helps in understanding behavior when the growing extracellular matrix has either 17 ruptured (due to hydrostatic forces) or it has created a cerebral edema in its vicinity. Hence if edema is 18 present, it would be better advised to consider the effect of fluid hydrostatic pressure. 19 Fourth, we have applied tumor forces on the walls of ventricles. This simplification allows us to 20 understand the biomechanical behavior and physiology of stenosis of CA. This approach is best suited 21 when considering intraventricular tumors. However, if the tumor is on any other location, then either the 22 interaction of brain parenchyma needs to be taken or a transformation matrix is needed to find the 23 transmission factor with which forces are either reduced or increased (as the case may be) as they travel through the brain parenchyma. Further, usually a force field or body forces may be enough to find the 1 major effect of tumor on walls of CA. However, sometimes stress field is also needed if brain 2 parenchyma interaction is considered. 3 Lastly, CSF biomechanical model is constructed by segmenting anatomically realistic brain ventricles. 4 We did not consider effect of Subarachnoid space (SAS) where 75% of CSF volume is usually found. The 5 purpose of not considering (SAS) in this study pertains to the fact that the biomechanical properties, flow 6 parametrization and Transmantle pressure distributions are largely driven or influenced by changes in 7 biomechanical behavior of CSF in ventricular regions apart from CA. For instance, obstructive 8 hydrocephalus is caused due to stenosis of CA thereby increasing ICP; changes of CSF velocity in 9 subarachnoid space are negligible because it is a large cavity, and thus creeping flow dominates. 10 Therefore, while analyzing CSF biomechanical properties, another dimension of studying the behavior of 11 CSF inside the ventricles may be a useful parameter. 12 13

CONCLUSION AND FUTURE WORK 14
This paper is to the best of our knowledge the first effort which model stenosis of CA under the 15 influence of brain tumor using FSI modelling. Our work shows that by using image informatics coupled 16 with FSI modeling, the underlying ventricular-CSF interaction can be well established and understood. 17 Furthermore, model proposed is not only restricted to the particular case at hand but can also be extended 18 to other cases where CSF-ventricular interaction is under consideration. For future work, we hope to 19 extend our proposed work and include particularly three aspects in our model. The factor of subarachnoid 20 space and brain parenchyma interaction which will define a complete intracranial environment during the 21 CSF-ventricular interaction. Second factor which can be included in future is the role of fluid forces of 22 tumor which may influence the outcome especially if they are accompanied by cerebral edema. Lastly, we 23 hope to incorporate in our model, role of blood-brain barrier and cerebral vasculature so as to obtain more The geometry is further refined by using Laplacian smoothing filter. Figure 1 shows axial scan of MRI 7 data and segmented ventricle and tumor body.

Meshing of Geometry and Mesh Independence Test 11
Mesh independence study is performed to find the ideal grid size. As Fluid-Structure Interface 12 technique is used therefore performing mesh independence test on fluid and structure domain separately 13 gives better indication of the converged solution. In solid domain, SHELL 181 elements were used. This simulation. Similarly, in the case of fluid domain, tetrahedral elements (SOLID 187) were used for 1 meshing. Mesh independence study was carried out to identify an ideal grid size. Results, confirmed that 2 a mesh density of 1144379 elements with minimum element size of 0.5 mm was good enough to capture 3 all the important information regarding peak velocity changes in the simulation. Figure 2 shows the 4 meshing of geometry and Fig 3. shows the mesh independence study for both solid and fluid domains.

Material Properties: Solid Domain 8
Choice of material properties for ventricular geometry was the most critical part of this study, 9 since it is responsible for overall model behavior. Density of ventricular geometry and the poisson ratio 10 was taken as 1000 kg/m 3 and 0.49 respectively [10]. The only thing that must be taken into knowledge 11 here is that, in literature a ventricular geometry is largely modeled as elastic when it is studied in an 12   to the availability of data from in-vivo studies. Prevailing data in the literature suggests that ventricles get 1 the properties of its surrounding tissues such as viscoelastic [13]. Since we are considering it in an 2 isolated environment, a linearized model was preferred in the literature [10]. Table 1 shows the material  3 properties. 4 5

Governing Equations 6
The numerical model on the fluid side consists of solving equations of conservation of mass and 7 momentum [14] given below: 8 (2) 10 Solving above equations (1) and (2) require spatial and temporal discretization together-with pressure-4 velocity coupling schemes. This can be used to interpolate pressures at the faces of a control volume. 5 Pressure-based solver is used with PISO (Pressure-Implicit with Splitting of Operators) coupling. For 6 spatial discretization, second order pressure, and for momentum second order upwind was used. For 7 gradient discretization, least square cell-based method was used. Implicit formulation is used in 8 discretizing equations. Implicit time marching scheme was employed using second order backward Euler 9 derivative formulation. 10 The governing equations relevant for the structural mechanics part are as follows: 11 where E is the stiffness matrix, is the stress vector and are the strains. The linear stress-strain law is 13 applied to get deformation on the body. Newmark Integration scheme is used to solve non-linear 14 equations and therefore implicit solver is needed to reach towards a solution. Implicit time marching 15 scheme is used together-with Newton-Raphson method for loads and displacement convergence. Newton-16 Raphson equation is given as follows 17 Equation (4) is used to converge forces, moments and displacements at each timesteps, which works by 19 the principle that the energy added due to the external loads must eventually balance the energy induced 20 by the reaction forces. Governing condition for FSI coupling is the kinematic coupling condition [15] in 21 which no slip condition with the fluid and structure interface was followed. The FSI equations for the 22 kinematic coupling at the interface Γ are given as follows (15): Here Γ ( ) is the displacement of the fluid mesh at the interface. The dynamic coupling is given by: 4 where is stress tensor on the fluid side n i is the normal vector and is stress tensor on the solid side. 6 The product . is called the traction vector. Due to displacements, dynamic meshing adaptability is 7 needed so that the moving cells on the fluid and solid side could be re-meshed correctly thereby 8 preserving shape quality. Diffusion based smoothing was used with diffusion parameter equal to 2. 9 10

Boundary Conditions and Tumor Growth Modeling 11
The boundary conditions were used to make the model as realistic as possible. Our model 12 consisted of two inlets and two outlets each of 3mm 2 diameter. These inlets were defined in the later 13 ventricles whereas outlets were defined in the fourth ventricle. We knew that the inside linings of walls of 14 ventricles have the Choroid Plexus; and in reality, the inside walls behave as a source for CSF. Whereas 15 in this study, inlets were defined at particular points (Fig 4), which is consistent with the previous studies 16 [7,10,12]. This is due to the fact that the behavior of CSF fluid flow is largely a creeping flow [7,10]. The 17 flow is consequently assumed to be in the laminar region throughout [7]. Therefore, in any case, defining 18 the entire inside walls of lateral ventricles as inlet has no ultimate bearing on velocity profile across the 19 CA. 20 1 2 At the fluid side, mass flow inlet boundary condition with bulk production of 500 ml/day (or 6.25x10 -6 3 kg/s) was used to mimic the real scenario of CSF flow [16]. While at the outlet, static pressure of 0 Pa 4 was defined [11]. CSF behaves very similar to water in terms of nature [15] with a density of about 1000 5 kg/m 3 and a constant viscosity of 0.00103 Pa.s. Whereas, the protein content has very limited and 6 negligible effect [15]. CSF was therefore, treated as Newtonian fluid. A Newtonian fluid is usually 7 described by Newton's law of viscosity presented in mathematical form as: 8 where is the shear stress, du/dy is the gradient of shear strain and is the dynamic viscosity constant. 10 Equation (9) suggests that if shear stress is directly proportional to rate of shear strain, then fluid will 11 behave as Newtonian fluid. This means that the shear stress along the walls would be maximum and the 12 fluid which is in contact with the medium has zero velocity/strain relative to the wall. This is hence, a no-13 slip condition. While the shear strain along the wall will be zero. 14 To add more value to this study, a novel schematic boundary condition was employed i.e. using CSF 1 pulsatility induced during cardiac cycle and its effect on the stenosis of CA. The pulsatile component is 2 ideally a linear combination of sinusoidal harmonics given as [17]: 3 Volumetric flow rate (t) = Bulk production + Asin(wt + α) + Bsin(2wt + β) (10) 4 In equation (10), bulk production, as noted above, is approximately equal to 0.3 ml/min (6.25x10 -6 5 kg/sec), A is equal to 0.21 ml/min (3.5x10 -6 ), and B is the part of second harmonic and is taken with an 6 amplitude of 0.05ml/min [17]. Carotid artery, which takes blood away from heart has zero phase  Regarding, the external loading which came from tumor, it is obvious that the initial tumor size needs to 13 be grown in temporal spectrum so that the loading conditions can become dynamic in nature. To predict 14 where F is the body force, ρ is the density of the tumor mass which is taken to be 1040 kg/m 3 [13], g is 11 acceleration due to gravity (9.81 m/s 2 ) and dV is the volume differential element. 1

Availability of data and materials 2
The datasets used and/or analysed during the current study are available from the corresponding author on 3 reasonable request. 4 We are thankful to HEC Pakistan, for providing us funding to execute this research. We are also indebted 1 to the National University of Sciences & Technology (NUST) for providing Computer Labs and logistic 2 support necessary to conduct this study. 3