Computational modelling of cerebrospinal fluid flow in a translational canine model

INTRODUCTION

Chiari malformation and syringomyelia are neurological disorders associated with cerebrospinal fluid (CSF) flow abnormalities which can lead to pain and neurological dysfunction. Studying CSF dynamics in affected human populations is limited by the low incidence and phenotypic variability of these conditions. In pediatric patients, type 1 Chiari malformation affects approximately 0.1% of the population, with concurrent syringomyelia present in about 30% of those cases^[1]. Population-level estimates are lower in Italy (0.008% for Chiari malformation and 0.005% for syringomyelia^[2]) and syringomyelia incidence in the Russian Federation is around 0.01%^[3], possibly reflecting regional familial predisposition.

In contrast, Chiari-like malformation, a naturally occurring analogue of human Chiari I malformation in dogs, and syringomyelia are highly prevalent and less phenotypically variable in the Cavalier King Charles Spaniel (CKCS). The incidence of Chiari-like malformation and syringomyelia in the CKCS ranges between 77%-100% and 15%-79%, depending on the population examined, exclusion criteria, etc. Although the statistical power of these studies is limited, these data show a significantly higher prevalence compared to humans, and increasing incidence with age^[4]. These high rates are linked to artificial selection for neotenous and brachycephalic cranial conformation, which reduces the volume of the cranial fossa and promotes neuroparenchymal overcrowding, especially at the cranio-cervical junction. Because these conditions arise spontaneously in the CKCS and share key pathophysiological features with the human disease, including posterior/caudal fossa crowding, altered craniospinal compliance, and syrinx formation, this canine breed has been proposed as a clinically relevant translational model for investigating disrupted CSF dynamics; providing access to a large volume of clinically acquired magnetic resonance imaging (MRI) data from affected CKCS which can be leveraged to study CSF flow. Our group alone has a dataset of hundreds of existing MRI sequences from which fluid simulations can be initiated. CKCS serve not only as experimental models, but also as clinical patients, such that mechanistic insights derived from this population have direct implications for veterinary diagnosis, surgical planning, and outcome assessment. To date, only 2 studies have modelled CSF flow in dogs^[5,6]; however, these primarily addressed mechanical stress in the spinal cord rather than the characterization of CSF flow patterns.

The use of dogs as translational models for CSF associated conditions is well established in comparable conditions. Early hydrocephalus experimental work relied heavily on canine models as their gyrification, ventricular anatomy, and size permitted direct surgical manipulation and physiological measurement that was comparable to humans. The first reproducible models of hydrocephalus are generally attributed to Walter Dandy and Kenneth Blackfan, who induced ventricular enlargement in dogs by obstructing CSF pathways and thereby demonstrated that hydrocephalus could be produced through interference with CSF circulation^[7-9].

Due to their gyrencephalic brains and size, dogs also served as key translational models for evaluating CSF diversion procedures, allowing implantation of shunt systems and assessment of long-term outcomes. Becker and Nulsen, 1968 used experimental dog models to evaluate valve-regulated ventriculo-venous shunts, demonstrating that pressure-controlled drainage improved long-term shunt function and reduced complications such as over-drainage and venous reflux^[10]. These principles were translated to veterinary patients, including ventriculo-atrial shunting in dogs^[11].

In addition to hydrocephalus, experimental work in dogs has been central to shaping modern concepts of syringomyelia as a disorder of CSF hydrodynamics and pressure transmission. One of the earliest reproducible large-animal models was described by Williams and Weller (1973), who produced syrinx-like cavities through intramedullary injection of CSF or saline in beagles. The development of these cavities supported a hypothesis through which fluid under pressure dissects through mechanically vulnerable regions of the spinal cord and may communicate secondarily with the central canal rather than arising from primary canal dilation^[12].

In contrast to syringomyelia and hydrocephalus, no equivalent experimentally derived canine translational models have been systematically developed for human Chiari malformation, however, there is sustained clinical interest in the translational relevance of naturally occurring Chiari-like malformation and syringomyelia in dogs. Particularly given the anatomical, pathophysiological, and imaging parallels with the comparable human conditions. Veterinary neurologists and researchers have regularly contributed to multidisciplinary clinical and scientific meetings such as the British Chiari Syringomyelia Group, reflecting recognition of the translational relevance of veterinary research^[13].

CSF flow in the craniospinal system is inherently pulsatile, exhibiting periodic variations in velocity and direction driven primarily by the cardiac cycle^[4]. In dogs, CSF circulates through a ventricular system that is continuous with the subarachnoid space (SAS), with caudal displacement during systole and rostral return during diastole. The anatomy of this system is shown in Figure 1. Several anatomical features distinguish the canine system from that of humans and are relevant to CSF hydrodynamics, including a patent central canal throughout the spinal cord, the absence of cerebellar tonsils, and the absence of the foramen of Magendie. The volume of CSF displaced into the spinal SAS depends in part on the compliance of the spinal dura mater^[14]; a measure of volumetric change per unit of applied pressure. These features collectively shape craniospinal coupling in dogs and must be accounted for when characterizing CSF flow dynamics. Additional CSF circulation is facilitated through the parenchyma via the glymphatic system and intramural peri-arterial drainage system; the examination of which is outside of the scope of this study. For a more comprehensive discussion of CSF anatomy, physiology, and flow in dogs, readers are referred to our recent review^[4].

Figure 1. (A) The cerebrospinal fluid (CSF) anatomy of the brain and spinal canal in the mid-sagittal plane; red indicates the choroid plexus, black arrows indicate flow directionality; (B) An example 3D model of the CSF anatomy of the brain and spine used in this study for computational modelling.

Phase contrast-MRI (PC-MRI) has been used in some studies to characterize CSF flow dynamics in dogs; a summary of these findings is provided in Table 1. Despite their widespread use, PC-MRI CSF velocity measurements are frequently treated as direct physiological surrogates without sufficient consideration of how anatomical location, flow complexity, and measurement methodology influence the reported values. This creates a critical translational challenge: variability in CSF metrics may reflect methodological artifacts rather than true disease-related differences, limiting their reliability as clinical or research biomarkers.

While PC-MRI provides valuable insights into CSF flow dynamics, it is inherently limited to in-plane velocity measurements at a small number of locations and at a relatively low spatial and temporal resolution. Out-of-plane flow components and direct pressure measurements cannot be captured with this technique, and therefore only a partial description of CSF dynamics is obtained. Novak et al. (1974) investigated CSF pressure in dogs using radioisotope cisternography and reported the CSF pressure waveform over multiple cardiac cycles with mean pressure 117 mmH₂O^[15]. Physiological mechanisms influencing CSF pressure dynamics, including transient pressure propagation associated with muscular stimulation, have also been demonstrated experimentally^[16], but quantitative data on the CSF dynamics in canines at rest over the cardiac cycle remain scarce.

Computational fluid dynamics (CFD) provides a complementary approach for resolving mechanistic uncertainties that limit the clinical translation of CSF flow measurements derived from in vivo imaging. CFD enables high-resolution quantification of velocity and pressure fields throughout the entire CSF space. Unlike PC-MRI, CFD allows interrogation of complex three-dimensional flow patterns, out-of-plane motion, and pressure gradients. This makes CFD particularly well suited for investigating how subtle anatomical variations influence global CSF dynamics. However, computational modelling studies of CSF dynamics in dogs are scarce and limited to semi-idealized models of the spinal column subject to hypothetical flow inputs^[5,6].

While CFD modelling of CSF in humans is a well-developed field^[17-21], comparable mechanistic characterization in the canine craniospinal axis is lacking, and it remains uncertain whether the key features of CSF flow identified in humans are conserved across species with distinct geometric scaling and distinct cranial morphometry. As a first step in addressing this gap, the present study provides the first CFD-based quantification of CSF hydrodynamics in the dog, simulating the entire cranio-spinal space and physiologically realistic boundary conditions consistent with arterially driven CSF flow; hypothesizing that CFD-based quantification of the CSF space can be achieved across multiple subjects using models derived from routine clinical MRI data. This study examines CKCS presenting without clinical signs of pain or other neurological dysfunction and scored Chiari-like malformation grade 0 or 1 on the Kennel Club British Veterinary Association Chiari-like malformation/syringomyelia (CMSM) health scheme^[22]. By establishing a validated baseline for CSF hydrodynamics in CKCS with no evidence of CSF associated abnormalities, this work begins to address a key methodological barrier to translation: the interpretation and standardization of CSF flow metrics used to guide clinical decision-making in veterinary patients. These initial unaffected models establish a methodological framework that can be readily extended to a cohort of CKCS with varying degrees of Chiari-like malformation/syringomyelia, leveraging the abundance of clinically acquired MRI data to systematically explore how progressive anatomical deformation influences CSF flow, how surgical methods can be optimized to improve long-term outcomes, and how drug delivery strategies can be refined by enabling tracking of solute transport, optimizing delivery location, infusion rate, and delivery mechanism from clinical imaging data.

MATERIALS AND METHODS

Boundary conditions and physiological assumptions

The simulations were designed to reproduce the oscillatory component of CSF motion driven by the cardiac cycle. Cardiac-driven expansion of the cerebral vasculature during systole transiently increases intracranial volume, displacing CSF caudally toward the spinal compartment, where compliance is greater. During diastole, this process reverses and CSF moves rostrally. This mechanism was represented by a time-varying velocity condition applied uniformly over the ventricular surfaces and the inner boundary of the cranial SAS. The imposed velocity waveform was based on experimentally measured CSF flow data in dogs and was constructed to satisfy key physiological constraints: pulsatility synchronized to heart rate, dominance of caudal flow during systole, temporal asymmetry between systolic and diastolic phases, and zero net flow over one cardiac cycle^[25,26]. This ensured that simulated CSF motion represented redistribution of fluid within the craniospinal system, rather than long-term accumulation or drainage.

To account for the compliance of the spinal compartment, the caudal end of the model incorporated a lumped representation of spinal dural compliance. This compliance allows temporary storage and release of CSF during the cardiac cycle and thereby modulating pressure gradients along the spine. This approach captured the functional effect of spinal compliance without explicitly modelling dural deformation. The boundary condition waveforms are provided in Figure 2.

Figure 2. Physiologically realistic velocity and pressure boundary conditions used in computational modelling of the CSF flow system over the cardiac cycle, detailed in [Supplementary Section 1]. CSF: Cerebrospinal fluid.

These assumptions isolate the oscillatory portion of CSF pressure while accurately modelling the dynamics of CSF movement.

RESULTS

Geometric and hydrodynamic characteristics of CSF flow in the canine craniospinal system were quantified. Flow simulations characterized pulsatile CSF dynamics, assessing regional flow rates, velocity distributions, Reynolds and Womersley numbers, and pressure gradients. An example velocity field outputted by a single model is shown in Figure 3.

Figure 3. Instantaneous velocity field for a single model at peak systole. The velocity is presented relative to the peak velocity in the spinal SAS, and truncated at 1. SAS: Subarachnoid space.

CSF velocity

Figure 4 shows the relative CSF a. max and b. mean velocity along the craniospinal axis normalized to the max/mean velocity and C2; a region commonly reported in CSF in vivo studies. Mean max velocities followed the same trend, with max relative velocities exhibiting greater amplitude. Relative velocities were lowest in the foramen magnum, before increasing to around 1 in the cervical SAS. A further increase in velocity was observed at the cervico-thoracic junction and continuing into the thoracic spine. In the central canal, as the anatomy was idealized to a pipe of approximately constant diameter, the variance in velocity was minimal down the spinal canal. Peak CSF velocity across the models was approximately 0.4 cm/s in this region.

Figure 4. (A) max and (B) mean velocity down the craniospinal axis for each CFD simulation, relative to the velocity at C2. Mean is presented with error bars of 1 standard deviation. CFD: Computational fluid dynamics.

Figure 5A illustrates the distribution of CSF velocity within the cerebral aqueduct, relative to the maximum velocity measured at a transverse aqueduct slice. Velocity contours revealed evidence of flow separation, jet-like structures, and secondary flow patterns within the mid‐ampulla of the aqueduct. Velocity profiles demonstrated marked deviations from parabolic or plug-like flow, with skewed and multi-peaked profiles indicative of complex flow behaviour [Figure 5B]. Figure 5Ci shows jet formation at peak systole, with flow recirculation around the dorsal aspect of the mid-ampulla of the cerebral aqueduct; in peak diastole, this region was observed to be near-stagnant [Figure 5Cii].

Figure 5. (A) Normalised CSF velocity magnitude within the cerebral aqueduct at peak systole (20% of the cardiac cycle), shown relative to the maximum velocity measured at a transverse slice through the mid-ampulla; (B) Flow rate through the transverse slice at the mid-ampulla of a single model shown as a red dashed line in (A), with the velocity profiles through the middle of the cross section shown at various points in the cardiac cycle; (C) Vector plot and streamlines of velocity magnitude in the mid-ampulla of the cerebral aqueduct at i. peak systole and ii. peak diastole relative to peak velocity in systole. CSF: Cerebrospinal fluid.

DISCUSSION

It is necessary to compare the results of a computational simulation to known quantities reported in the literature to assess the physiological plausibility of the results. To date, no studies have directly measured CSF pressure waves in dogs; hence, this study is limited to CSF flow velocities reported across a range of breeds. Such studies generally report velocities from broader anatomical regions (e.g. “foramen magnum” or “cisterna magna”) and across mixed breeds and body sizes. Where possible, data from the literature were categorized based on the anatomical descriptions and figures provided, to allow comparison at corresponding locations in the models.

Christen et al. (2021) and Cho et al. (2021) investigated CSF flow within the cerebral aqueduct and cervical spinal SAS in cohorts of six Beagles^[27,28]. Both studies reported pulsatile, bidirectional CSF flow with marked waveform asymmetry, characterised by narrow, high-amplitude systolic peaks and broader, lower-amplitude diastolic troughs. A ventral-dominant flow distribution was also observed, with reduced dorsal motion, likely attributable to SAS geometry and posture during imaging. Cerda-Gonzalez et al. (2009) examined CSF flow in the spinal SAS of fifty-nine CKCS and reported more complex flow patterns^[29]. Localised regions of high velocity were identified and were strongly correlated with syrinx presence and size, although not with neurological clinical sign severity. Dorsal CSF velocities were significantly reduced in CKCS; however, the degree of ventral flow dominance was lower than that reported in Beagles^[27,28], likely reflecting morphological differences between brachycephalic CKCS and mesocephalic Beagles. Rich et al. (2024) investigated physiological variability in a population of 32 neurologically normal dogs of mixed breeds, of which four were classified as brachycephalic; finding that body weight, age, and sex significantly influence peak CSF velocities^[25].

It should be noted that the aqueduct velocities reported by Farke et al. (2024) were considerably higher than those in the aforementioned studies^[30]. This discrepancy is likely due to the more rostral sampling location selected at the junction between the cerebral aqueduct and the third ventricle. At this level, the aqueduct narrows substantially, and, due to conservation of volumetric flow, measured velocities increase accordingly.

Z-scores were calculated using the reported literature means [Table 1] and standard deviations, providing a measure of deviation relative to the observed inter-subject variability in each study. Given the anatomical and methodological variability inherent in PC-MRI CSF velocity measurements, z-scores within approximately ±2 were considered to indicate reasonable agreement. These data are presented in Table 4.

Using the weight-stratified datasets reported by Rich et al. (2024), CFD velocities at the insertion site of the apical ligament yielded z-scores of -1.2 for dogs < 10kg, -0.8 for dogs between 10 kg and 20 kg, and -1.0 for the overall dataset^[25]. At C2, corresponding z-scores were -0.5, -0.9, and -1.1, respectively. These values indicate that the CFD-derived velocities lie well within 1 to 1.5 standard deviations of the reported experimental means, suggesting good agreement given the differences in breed and methodology.

Comparisons with data from beagles revealed larger deviations at the insertion site of the apical ligament, with z-scores of -3.4 relative to Cho et al. (2021) and -1.3 relative to Christen et al. (2021)^[27,28]. The moderately large deviation from the data presented by Cho et al. (2021) reflects a combination of breed-specific anatomical differences, small sample size, and the small relative size of their standard deviation^[28]. In contrast, the comparison with Christen et al. (2021) yields a z-score in closer agreement. At the dens-axis junction, comparison with Christen et al. (2021) resulted in a z-score of -1.8, again consistent with moderate deviation within the expected physiological spread^[27].

Using data presented by Cerda-Gonzalez et al. (2009), CFD-derived peak velocities produced z-scores of -0.03 relative to controls and 0.5 relative to unaffected CKCS at the dens-axis junction, indicating very good agreement with both cohorts. At the C2/3 level, however, comparison with unaffected CKCS resulted in a z-score of 6.9. This large deviation appears to be driven primarily by the small reported standard deviation in Cerda-Gonzalez et al. (2009) combined with systematically lower peak velocities than those reported in most other studies. While breed-specific differences may contribute, this is unlikely to fully explain the discrepancy, as the present simulations were also based on CKCS anatomy. Methodological factors may therefore play a role; Cerda-Gonzalez et al. (2009) employed earlier phase-contrast MRI techniques with lower spatial and temporal resolution than more recent studies, which may alias peak velocity estimation^[29].

These data demonstrate that substantial variability exists in reported CSF velocities within the literature, particularly for measurements broadly attributed to the foramen magnum or cisterna magna. Importantly, large z-scores do not necessarily indicate poor agreement but rather reflect the sensitivity of peak velocity measurements to precise anatomical location, sampling strategy, and cohort characteristics. The alignment with literature values supports the physiological plausibility of the simulated flow fields and provides confidence in using the model to explore spatial flow patterns along the spinal canal.

The results indicate that within the cisterna magna alone, peak systolic velocity varies markedly between closely spaced anatomical landmarks. Consequently, studies that report CSF velocity measurements from the cisterna magna or foramen magnum without precise anatomical specification may inadvertently conflate distinct flow regimes. These findings highlight the importance of rigorous anatomical definition when acquiring and interpreting CSF velocity data, particularly when comparing results across studies or using such measurements for model validation.

Another commonly reported location for CSF velocity measurement in veterinary studies is the cerebral aqueduct; most frequently sampled at the mid-ampulla^[25,27,28], or at the rostral end of the aqueduct connecting to the third ventricle^[30]. These locations are often selected due to their clear anatomical definition.

In the present models, the cerebral aqueduct exhibited markedly complex flow behavior. As shown in Figure 5, peak systolic flow was characterized by flow separation, jet-like structures, and pronounced secondary flow patterns, particularly within the distended mid-ampullary region. Velocity profiles demonstrated substantial spatial heterogeneity, with localized high-velocity jets adjacent to regions of recirculation indicating that the cerebral aqueduct does not support a simple, fully developed flow regime.

These complex flow patterns are likely generated by abrupt changes in cross-sectional area along the aqueduct, especially the transition between the relatively constricted rostral segment and the expanded mid-ampulla. In brachycephalic breeds, and particularly in dogs with neuroparenchymal overcrowding, constriction of the rostral aqueduct may be exacerbated, increasing the likelihood of jet formation, flow separation, and secondary flow patterns. Under these conditions, local peak velocities become highly sensitive to both anatomy and sampling location.

These findings suggest that the cerebral aqueduct represents an inherently unreliable site for quantitative CSF velocity sampling. While anatomically convenient, the highly complex and spatially variable flow field makes reported velocities strongly dependent on precise measurement location and temporal resolution. Measurement techniques sampling below the time-/length-scales upon which complex flow patterns occur are therefore particularly vulnerable to spatial and temporal aliasing in this region. As a result, cerebral aqueduct velocity measurements should be interpreted with caution, especially when used for cross-study or cross-cohort comparisons.

Table 3 shows the regional partitioning of CSF flow, with a significantly higher portion of CSF volume exchange occurring in the SAS, consistent with its larger cross-sectional area and lower hydrodynamic resistance. This implies that, in the first approximation, central canal is negligible and can therefore be neglected. However, syrinx presentation in dogs is predominantly communicating, with a patent central canal. Therefore, from a syringomyelia research perspective, flow in the central canal is physiologically significant and therefore should be represented in canine CSF models.

The Womersley numbers in the SAS and central canal indicate that CSF flow exhibits a range of pulsatile behavior along the spinal axis. In the SAS, α ranged from 2.0 to 7.6, suggesting that inertial effects become more pronounced likely reflect cross-sectional tapering leading to stronger transient inertia effects. In contrast, the lower Womersley numbers in the central canal (< 0.5) indicate that viscous effects dominate, and velocity profiles are broadly parabolic throughout the cardiac cycle.

Velocity profiles in the SAS displayed some characteristic behavior for pulsatile inertia-dominated flow, with a central “plug” that lagged behind near-wall velocities, whereas flow in the central canal was more uniform and dominated by viscous diffusion. This distinction has implications for solute transport and mixing, as areas with higher α may experience greater phase lag between pressure and velocity, enhancing shear-driven dispersion. These findings were consistent with the low Reynolds numbers observed throughout the domain (< 50), confirming that CSF motion remained laminar and strongly influenced by viscous forces. Together, the Womersley and Reynolds number distributions suggest that while unsteady inertial effects shape the pulsatile response of the SAS; the overall regime in both compartments remained viscous-dominated. Compared with previous studies in humans and monkeys, the Womersley numbers observed in the dog are lower, further indicating a greater relative influence of viscosity in CKCS CSF flow^[21,31-36]. It is worth noting that CSF in the CKCS, a brachycephalic breed, will exhibit different flow patterns to a breed with a comparatively larger SAS.

Figure 4 demonstrates that peak CSF velocities reached values up to approximately three times those measured at C2, exceeding the magnitudes reported in previous canine MRI studies. However, existing experimental studies typically report CSF velocities caudally to the C2/3 intervertebral disc, and CSF flow characteristics beyond this level have not been described. The elevated velocities observed in the present models therefore do not represent disagreement with the published literature but rather reflect an extension of CSF flow characterisation into more caudal regions of the canine spine that have not previously been quantified experimentally.

Figure 4 also reveals spatial trends in SAS velocity along the craniospinal axis. Rostral to C2, CSF velocities are relatively low, increasing towards the cervical spine and remaining broadly consistent between C2 and C7. A further increase in velocity is observed at the transition from the cervical to the thoracic spine.

When maximum velocity was extracted independently at each transverse plane, a consistent alternating pattern was observed along the cervical spine, with lower peak velocities at intervertebral disc levels and higher peaks at the midpoints of vertebral bodies. Because the flow is pulsatile and peak values correspond to nodal maxima, the observed variation is plausibly explained by spatial variations in local hydraulic impedance within the SAS cross-section, which introduce phase differences in the timing and location of peak velocity. This effect is substantially attenuated in the mean velocity profiles, supporting the interpretation that the observed modulation in peak velocity arises from localised, in-plane impedance variations.

In the present models, a transmural pressure gradient was observed between the central canal and SAS across the spinal cord. This was investigated further using a single model representing the spine from the foramen magnum to the level of the second lumbar vertebra. This indicated that, while these results indicate that such a gradient can arise under physiologically realistic pulsatile forcing, the nature of this gradient was found to be highly sensitive to the hydraulic impedance of spinal segments not represented in the truncated domain. Both the presence and character of the pressure gradient depend strongly on the downstream impedance. Consistent with this sensitivity, the central canal and SAS exhibited distinct temporal pressure and flow responses to the same pulsatile boundary conditions. The central canal’s high viscous resistance produced a rapid response to pressure reversal, followed by a delayed secondary adjustment as fluid volume was redistributed between the two pathways. Future investigations should specifically target this behavior using models explicitly designed to resolve full-length spinal impedance and inter-compartmental coupling, as transmural pressure gradients may play a critical role in the pathogenesis of syringomyelia.