In general, we consider a perivascular tree-like domain Ω consisting of a network of branching generalized annular cylinders Ω ⁱ , with Ω ⊆∪ _i∈I Ω ⁱ , spatial coordinates x ∈ Ω and time t ≥ 0. The boundary is denoted ∂Ω, with boundary normal n. We assume that each generalized annular cylinder Ω ⁱ has a well-defined and oriented, topologically one-dimensional centerline Λ ⁱ with coordinate s. We set Λ = ∪ _i∈I Λ ⁱ . Along s, we define the cross-sections C ⁱ = C ⁱ (s, t) of Λ ⁱ with area A ⁱ = A ⁱ (s, t). We denote the inner radius of Ω ⁱ by R 1 i and the outer radius of Ω ⁱ by R 2 i ; these radii will in practice vary with s, t and the angular coordinate θ. We denote the set of bifurcation points i.e. the points at which the centerlines of branches meet by B .

We introduce three specific geometries of increasing complexity: from an axisymmetric cylinder (A) to an image-based perivascular geometry without any bifurcations (B) and one with a bifurcation (C) (Figure 1 and Table 1). The image-based geometries (B) and (C) are constructed from non-pathological artery segments from the Aneurisk dataset repository [30], and thus define high-fidelity 3D representations of human brain surface arteries. In each of these geometries, the PVS domain is defined by creating a generalized annular cylinder surrounding the vascular segment with the vascular wall as the inner surface of the PVS. The width of the PVS is set proportional to the blood vessel diameter (by factor of 0.95) and scaled (to a mouse scale) [18, 31]. Three-dimensional PVS flow in geometries A and C have been studied previously [18] and will be used for comparison. We define as PVS inlets and outlets (∂Ω_in and ∂Ω_out) the PVS ends surrounding the vascular inlets and outlets, respectively, noting however that fluid may flow both in and out of both the inlet and outlets. We denote the inner PVS wall (boundary) by ∂Ω_inner and outer wall by ∂Ω_outer.

The 3D PVS construction and the 1D centerline extraction are performed using PVS-meshing-tools [32], largely based on VMTK [33]. The extracted centerline comes with underlying data including the branch lengths and vessel radii. The centerline radius refers to the radius of the maximal inscribed circle of the vessel cross-sections. The meshing of both 3D and 1D PVS domains is performed within PVS-meshing-tools [32] using meshio [34] and GMSH [35]. The centerline meshes consist of topologically one-dimensional intervals embedded in three dimensions. The bifurcation points b ∈ B ⊂ Ω are explicitly labeled within each centerline mesh. Each branch is also separately tagged and given a consistent orientation. This procedure allows for the identification of bifurcation points as the outlet of one (parent) centerline and the inlet of other (daughter) centerlines, and a split of the full perivascular network into oriented mesh branches.

Flow of CSF in surface PVSs is reported to be laminar, with low Reynolds numbers (10^–4−10^–2) and moderate Péclet numbers (10²–10⁴) for 1 μm spherical particles transported at low Reynolds number [31], a mean flow speed of up to 60 μm/s, and parabolic flow profiles [31]. We therefore model the flow of an incompressible, viscous fluid flowing at low Reynolds and Womersley numbers via the time-dependent Stokes equations over a time-dependent domain Ω = Ω(t) representing the PVS. The fluid velocity v = v (x, t) for x ∈ Ω(t) at time t and the CSF pressure p = p (x, t) then solve the following system of time-dependent partial differential equations (PDEs) [18, 36]: ρ ∂ t v − μ ∇ 2 v + ∇ p = 0  in  Ω t , (1a) ∇ ⋅ v = 0  in  Ω t , (1b) where ρ is the fluid density and μ is the dynamic fluid viscosity. To model CSF at body temperature, we set the fluid density to ρ = 10³ kg/m³ and the dynamic viscosity to μ = 0.697 × 10^–3 Pa s. As in our previous full models of perivascular flow [18], the initial PVS mesh defines the reference domain Ω(0), and we assume that Ω(t) at time t > 0 is given by a deformation d of the reference domain: Ω(0)↦Ω(t) with x = d (X, t), X ∈ Ω(0), x ∈ Ω(t). We denote the domain velocity associated with d by w (thus d ̇ = w ).

At the PVS ends, we prescribe a traction condition corresponding to a known, applied pressure p ̃ = p ̃ ( x , t ) : σ n ≡ μ ∇ u − p I ⋅ n = − p ̃ n  on  ∂ Ω in  and  ∂ Ω out . (2)

We either prescribe (i) zero pressure at both ends p ̃ = 0 , or (ii) a constant-in-time pressure gradient Δ p ̃ > 0 by setting p ̃ in = L out Δ p ̃ at the inlet, letting p ̃ out = 0 at the outlet furthest from the inlet with distance L _out, and setting p ̃ out at any other outlets such that the average pressure gradient over each branch path ( p ̃ i − p ̃ out ) / L out is constant and equal to the prescribed pressure gradient Δ p ̃ mmHg/m. This static pressure difference can represent e.g., a hydrostatic pressure difference, a venous pressure differential, or some other systemic pressure difference. Other types of boundary conditions could also be considered, see e.g., a discussion of compliance conditions in [18], or [37].

On the inner and outer PVS walls (along the length of the PVS), we set the fluid velocity v to match a known, prescribed domain velocity w = w (x, t). For the inner PVS wall, we either (i) consider a rigid wall and set v = w = 0, or (ii) impose a pulsating wall displacement: d | ∂ Ω inner X , t = A X , t n , (3) with reference to the initial (fixed) mesh with coordinates X and prescribe v = w = d ̇ . To represent wall motion induced by the cardiac pulse wave, we let the amplitude A be defined by the juxtaposition of an experimentally-observed wall motion time series [31] either applied uniformly along the length of the PVS or as a travelling wave along the PVS length with wave speed c = 1 m/s and frequency 10 Hz. We refer to [18] for the detailed description. To represent wall motion due to vasomotion, we consider a similar set-up but with a travelling sinusoidal wave in time with a frequency of 0.1 Hz and wave length λ = 8 mm [38], and an amplitude A of 7.5% of the initial inner radius R ₁. We note that for all models, the wall moves in the normal (radial) direction only. For the outer PVS wall ∂Ω_outer, we set v = w = 0.

The system starts at rest with v = w = 0 at t = 0. The system reaches the periodic steady state nearly immediately, and we report results starting from the first cycle.

We define a reduced, topologically one-dimensional, model approximation of the full PVS flow model [(1) with the given boundary and initial conditions] under the following stipulations. For each branch Ω ⁱ (t) with centerline Λ ⁱ and local coordinate system (s, r, θ), where s represents the path length (or axial coordinate), r is the radial coordinate and θ is the angular coordinate, we suppose that:

(I) Axial symmetry. Fields and input parameters are independent of the angular coordinate θ;

For the velocity profile v _vp, we here choose a normalized annular Poiseuille flow: v vp r = v poise r v poise R 1 + R 2 2 , v poise r = 1 − r 2 R 1 2 + R 2 2 − R 1 2 R 1 2 ⁡ ln R 2 / R 1 ln r / R 1 . (5)

This velocity profile is parabolic in r (as for Poiseuille flow in a cylinder) with a logarithmic correction that accounts for the annulus.

In particular, the domain velocity w is assumed independent of the angular coordinate θ. Note that we do not assume other velocity components (than the axial) to necessarily be zero. We emphasize that these assumptions will in general not be satisfied by realistic geometries and flows. Thus, the reduced model defines a model approximation associated with a certain modelling error.

Under the assumptions (I–V), the full PVS flow model can be reduced to the following system of time-dependent differential equations: find the cross-section flux q ̂ = q ̂ ( s , t ) and the cross-section average pressure p ̂ = p ̂ ( s , t ) such that for each centerline Λ ⁱ (denoting q ̂ | Λ i = q ̂ i and p ̂ | Λ i = p ̂ i ): ρ A i ∂ t q ̂ i − μ A i ∂ s s q ̂ i + μ α i A i q ̂ i + ∂ s p ̂ i = 0  on  Λ i , (6a) ∂ s q ̂ i = f ̂ i  on  Λ i , (6b) hold. f ̂ i s ≡ 2 π R 1 i s , t w R 1 , s , t ⋅ n | ∂ Ω inner + 2 π R 2 i s , t w R 2 , s , t ⋅ n | ∂ Ω outer . (7)

Moreover, A ⁱ = A ⁱ (s, t) denotes the cross-section area, while α ̂ i = α ̂ i ( s , t ) is a lumped flow parameter that depends on the domain geometry and the choice of velocity profile v _vp: α i s , t ≡ 1 A v ̄ ̄ vp s 2 π R 1 i s , t ∂ r v vp R 1 i s , t − 2 π R 2 i s , t ∂ r v vp R 2 i s , t , (8) and where v ̄ ̄ vp is the velocity profile integrated over each cross-section: v ̄ ̄ vp ≡ ∫ C s v vp r d r d θ . (9)

We also define the (one-dimensional) normal stress induced by q ̂ and p ̂ : σ ̂ ≡ μ A ∂ s q ̂ − p ̂ , (10) which corresponds to an average of the axial (s-)component of the normal stress in (2) over each cross-section

At the bifurcation points b ∈ B ⊂ Ω , we impose the following two conditions representing conservation of flux and continuity of normal stress, respectively: q ̂ p s p = q ̂ d 1 s d 1 + q ̂ d 2 s d 2 , (11) σ ̂ p s p = σ ̂ d 1 s d 1 = σ ̂ d 2 s d 2 , (12) where Λ ^p and Λ d 1 , Λ d 2 represent the centerlines of the parent and two daughter branches, respectively, associated with the bifurcation point b and s ^⋅ = ι ^⋅(b) where ι ^⋅ denotes the map from three-dimensional bifurcation point to the one-dimensional centerline coordinate for each branch Ω^⋅.

System (6) defines a set of equations for each branch centerline Λ _i and is closed by the bifurcation conditions (11, 12), together with boundary conditions at the PVS inlet and outlets, as well as initial conditions for the cross-section flux. Specifically, in place of the traction condition (2), we prescribe the corresponding pressure difference for the (average) normal stress σ ̂ cf. (10). In this manner, the (one-dimensional) solutions q ̂ i and p ̂ i of the reduced model (6) define approximations of the (three-dimensional) axial flux and pressure solving (1) integrated or averaged over each cross-section: q ̂ i s ≈ ∫ C i s v s s , r , t r d r d θ ≡ A i s q ̄ i s , p ̂ i s ≈ 1 A i s ∫ C i s p s , r , t r d r d θ .

The factor r originates from integrating in cylindrical coordinates. We note that the wall velocity w, which defines a boundary condition for the full PVS model (1), enters as a body force in the reduced model (6).

We solve the full PVS (1) via a previously developed and verified arbitrary Lagrangian-Eulerian (ALE) formulation and finite element discretization [18]. This solver builds on the standard FEniCS finite element software suite [39], and is openly available [40].

To compute numerical solutions to the reduced model (6), we consider a first-order implicit Euler scheme in time and a higher-order finite element method in space. The finite element mesh T of the centerline Λ is composed of mesh segments T i , one for each centerline branch Λ ⁱ . Each mesh segment is a mesh consisting of intervals embedded in R 3 . We label the set of bifurcation points B , inlet points I and outlet points O , and define the following finite element spaces:

• The flux space V _h is the space of continuous piecewise quadratics over T i for each i.

The flux is thus solved on each mesh segment representing the PVS network branches and may be discontinuous across bifurcations. We impose the flux conservation condition (11) weakly using a Lagrange multiplier formulation. The pressure is solved on the whole mesh and is continuous at bifurcations by construction.

For each discrete time t ^k , given q ̂ h k − 1 at the previous time t ^k−1 and time step Δt = t ^k − t ^k−1, we solve for the approximate cross-section flux q ̂ h k ∈ V h , average pressure p ̂ h k ∈ Q h and a Lagrange multiplier [corresponding to the normal stress (10) at the bifurcation points] λ h k ∈ R h solving a q ̂ h k , p ̂ h k , λ h k , ψ , ϕ , ξ = L k ψ , ϕ , ξ , (13) for all finite element test functions ψ ∈ V _h , ϕ ∈ Q _h , and ξ ∈ R _h . The left-hand side bilinear form a is defined by: a q , p , λ , ψ , ϕ , ξ = ∑ i ∈ I ∫ Λ i 1 A i ρ + Δ t μ α i q i ψ i + Δ t μ A i ∂ s q i ∂ s ψ i + ∂ s q i ϕ i − Δ t ∂ s ψ i p i d s + ∑ b ∈ B λ b ψ b + ξ b q b , (14) where λ ^b (or ξ ^b ) is simply the entry of the vector λ (or ξ) corresponding to bifurcation point b, and we define the natural jump: ψ b = ψ p b − ψ d 1 b − ψ d 2 b . (15)

The right-hand side linear form L is: L k ψ , ϕ , ξ = ∑ i ∈ I ∫ Λ i ρ A i q ̂ h k − 1 , i ψ i + f i ϕ i d s − ∑ x ∈ I Δ t p ̃ in x ψ i I x + ∑ x ∈ O Δ t p ̃ out x ψ i O x , (16) where the superscript i _I (i _O ) in the inlet (outlet) terms above refers to the unique centerline branch associated with the inlet (outlet) points.

The numerical solver for the reduced model was implemented in the well-established FEniCS Project finite element software [39]. The solver, and in particular the definition of the partially continuous flux space, builds on mixed-domain features [41] and relies on the latest development version of FEniCS. All data and source code are available via Zenodo [42].

An overview of the six computational models considered is given in Table 2. Each model is labeled with reference to its domain (A, B, or C) followed by a number indicating the driving forces included: (1) a given pressure drop, (2) wall movement due to cardiac pulsations and (3) wall movement due to vasomotion. For each model, we consider the full three-dimensional version as well as the reduced model.

To compare the solutions from the full and reduced models, we consider the following quantities of interest. For each domain, we define a set of cross-sections as follows. For domain A, we define the left-most end as the inlet (s = 0) and define an upper cross-section. For domain B, we consider the inlet and outlet ends of the PVS, as well as upper and lower cross-sections. For domain C, we consider the inlet at s = 0, and the two outlets, as well as three additional cross sections near the inlet, on the largest daughter branch relatively close to the bifurcation, and near the outlet of the other daughter branch.

We then compute for each cross section C(s) the averaged pressure p ̄ ̄ h ( s , t ) and the cross-section flux q ̄ ̄ h ( s , t ) : p ̄ ̄ h s , t = 1 | C | ∑ k w k p h x k , t , q ̄ ̄ h s , t = ∑ k w k v h x k , t ⋅ n C x k (17) for a quadrature scheme with points x _k and weights w _k defined over C and an approximation |C| of the cross-section area. Here n _C is the normal vector of the cross-section. The averaging is implemented by using the Frenet frame associated with Λ to map from an annular cylinder in a reference domain onto the cross-section, similar to the implementation of the averaging operator in fenics_ii [43].

With this in hand we define the percentagewise relative model discrepancy E _q (t) in the flux by E q t = ‖ q ̄ ̄ h t − q ̂ h t ‖ L 2 Λ ‖ q ̄ ̄ h t ‖ L 2 Λ ⋅ 100 % (18) and similarly for the pressure E _p . We typically compute this quantity if the flow is driven by a constant pressure gradient. In this case the fluid starts at rests and then quickly develops to stationary, annular Poiseuille flow. We then compute E _q (T) and E _p (T), where T denotes the final time.

For pulsatile flow, we typically compare the percentagewise relative error in peak pressure e p peak ( s ′ ) and peak cross-section flow e q peak ( s ′ ) at some cross section C (s′, t), where e q peak s ′ = | max t ∈ 0 , T q ̂ h s ′ , t − max t ∈ 0 , T q ̄ ̄ h s ′ , t | | max t ∈ 0 , T q ̄ ̄ h s ′ , t | ⋅ 100 % (19) and e p peak ( s ′ ) is similarly defined. Finally, we compare the net fluxes Q of the full and reduced model, where Q associated with the velocity v = v (x, t) can be computed as: Q = ∫ 0 T ∫ ∂ Ω in v ⋅ n d x d t , (20) and the corresponding quantity associated with the flux q ̂ = q ̂ ( s ) by: Q = ∫ 0 T ∑ x ∈ I q ̂ x d t (21) where the integration in time is over one period [0, T].

Flow in an axisymmetric annular cylinder of length ℓ driven by a constant pressure difference Δp (Model A1) is described by the analytic expression: q ̂ s , t = A Δ p μ α ℓ 1 − exp − μ α t ρ , p ̂ s , t = Δ p ℓ s + p ̂ 0 , (22) where α is the lumped flow parameter given by (8) and which is constant in time and space in this case. For the velocity profile (5) defined over geometry A (cf. Table 1), α = 7325.3/m², and μα/ρ = 5105.7/s. Thus, the time-dependency is negligible after only a few milliseconds, and the flow develops near-instantaneously to steady-state Poiseuille flow.

Both the full and reduced models reproduce the exact annular Poiseuille flow characteristics of this case (Figure 2A). The numerical difference between the analytic and computed reduced solutions for the cross-section flux q ̂ and average pressure p ̂ is at machine precision ( ‖ q ̂ ( T ) − q ̂ h ( T ) ‖ = 1.7 × 1 0 − 14 and ‖ p ̂ ( T ) − p ̂ h ( T ) ‖ = 2.6 × 1 0 − 17 ) (T = 1 s). In general, the total error is the sum of the model error and the numerical error associated with the space-time discrete approximation (13). For Model A1, the model error is zero as the model reduction assumptions (I–V) are exactly fulfilled by the geometry and flow pattern. As the total error also vanishes, we note that the numerical error is also negligible for this case.

Next, we examine the PVS flow and pressure generated by uniform axisymmetric pulsations of the inner PVS wall (Model A2, Figures 2B–E). The inner wall movement changes the inner domain radius R ₁ in time. The fluid is pushed out at the both ends as the PVS width decreases, and flows back in at both ends as the PVS width returns to baseline. This behaviour is reproduced by both the full (Figure 2B [18]) and reduced models (Figure 2C). We note that the reduced model assumptions (IV, V) do not hold in this scenario as the PVS axial velocity profile is no longer identical to the Poiseuille velocity profile, and the pressure is not perfectly constant on each cross-section. Comparing the full and reduced cross-section fluxes q ̄ ̄ h and q ̂ h , we observe however that the two models still agree closely (Figures 2D,E), both at the inlet and at an interior cross-section. Moreover, the time-profile of the reduced and full cross-section flux approximations are very similar (both at the inlet and at the interior cross-section, Table), though with small (Δt s) shifts in time. The peak outfluxes occur at the inlet and outlet; the peak outflux for the full model is 1.53 × 10^–3 μm³/s, and 1.47 × 10^–3 μm³/s for the reduced model (Figure 2D). The peak pressure occurs in the middle of the domain; the peak pressure for the full model is 0.194 and 0.193 Pa for the reduced model. Using (19) the relative error in the peak cross-section flux at the inlet is e q peak ( 0 ) = 4.0 % and in the peak cross-section (average) pressure e p peak ( 0 ) = 0.8 % . There is thus a small discrepancy between the two models, as expected by the violation of the reduced model assumptions.

In contrast to the axisymmetric geometry A, the image-based geometries B and C express angular and axial variations in radius. The inner and outer radii of these geometries vary along the length of the domain (with s) and depend on the angular coordinate θ, with the latter violating model assumption I. To study the resulting model error in isolation, we again examine the pressure-driven flow predicted in full and reduced models but now of geometry B (Model B1, Figure 3). The full numerical approximation of the pressure is nearly constant over each cross-section. On the other hand, the velocity profile varies between cross-sections and with the angular coordinate within each cross-section (Figure 3A). Therefore, we expect a larger model error in the reduced model compared to the previous case(s). At steady state (t = 0.5), the reduced pressure approximation p ̂ varies nearly linearly along the length of the domain as expected, and the reduced flux approximation q ̂ is essentially constant along the centerline with value q ̂ = 4.28 × 1 0 − 4 μ L/s. Computing the corresponding cross-section flux from the full model, we find values ranging from 3.5 × 10^–4 to 5.31 × 10^–4 μL/s. The relative model discrepancy (18) in the pressure is E _p = 2.6% and for the flux E _q = 12.6%.

Cardiac wall motion and vasomotion may drive pulsatile perivascular flow with different flow characteristics. To evaluate the model discrepancy induced by different physiological drivers, we compare the full and reduced models over an image-based PVS segment driven by wall motion induced by the cardiac pulse wave (Model B2) and by vasomotion (Model B3). The cardiac pulse wave induces wall motion at a higher frequency (10 Hz) travelling at a higher wave speed (1000 mm/s), while vasomotion creates pulsations at lower frequencies (0.1 Hz) and at a lower wave speed (0.8 mm/s). Both models include angularly, axially and temporally varying radii, and we expect model assumptions I, IV-V to not hold.

Both pairs of models induce pulsatile bidirectional flow in and out of the PVS segment in synchrony with the pulsating wall (Figure 4, Supplementary Video S1) with peak pressure magnitude in the middle of the segment, and conversely, low velocities in the middle of the domain and higher velocities near the PVS ends. Both model scenarios lead to pressure fields that are nearly constant on each cross-section (Figures 4C, 5), but with angularly varying velocity profiles (Figures 4B, 5).

For the cardiac wall motion, the overall cross-section average of the full pressure p ̄ ̄ h ranges from −0.05 to 0.26 Pa, while the full cross-section flux v ̄ ̄ h ranges from −1.54 × 10^–3 to 1.95 × 10^–3 μL/s. The reduced model accurately captures the temporal and spatial characteristics of the full model (Figures 4D–G). For the reduced model, the overall cross-section pressure p ̂ h ranges from −0.06 to 0.29 Pa, while the cross-section flux q ̂ h is between −1.61 × 10^–3 and 2.23 × 10^–3 μL/s. Comparing the full and reduced pressure and flux over time at an interior, lower cross-section with axial coordinate s′ (Figures 4B,C), we observe that the reduced model slightly overestimates the peak pressure and flux when compared to the full model (Figures 4F,G). Using (19) at this cross-section we find that the relative error of the peak pressure is e p peak ( s ′ ) = 19.0 % and the relative error of the peak flux is e q peak ( s ′ ) = 1.2 % . One shall note that the space discretization of the initial 3D model has a non negligible impact on these relative errors. Indeed, using a finer 3D mesh composed of 333000 tetrahedrons instead of the initial 63000 lowers the relative error of the peak pressure e p peak ( s ′ ) to 12.1%. The relative error on peak flux is not significantly impacted but e q peak ( s ′ ) was already very small.

For the vasomotion scenario, the domain movement is larger compared to the cardiac wall motion, but the wall velocity is lower (peak wall speed of 0.001 vs. 0.005 mm/s). The resulting peak (in terms of magnitude) cross-section pressure is −0.012 Pa and peak cross-section flux is 9.14 × 10^–5 μL/s (Figure 5, Supplementary Video S2). These are one-to-two orders of magnitude lower than for the cardiac wall motion scenario. Comparing the full and reduced models in two interior (upper and lower) cross-sections, we observe that the cross-section pressure q ̂ h matches pulsatile behaviour of the average cross-section pressure in the full model q ̄ ̄ h (Figure 5B) but that the peak amplitude is higher. The largest model differences in pressure at the lower cross-section is at the peak pressure; there the relative difference in peak pressure is e p peak ( s ′ ) = 52.7 % . The similar observations hold for the flux, but the model discrepancies are lower: the relative difference in peak flux is e q peak ( s ′ ) = 15.6 % . Moreover, the full and reduced models agree on a pressure phase shift of 0.5s. In agreement with our previous findings, the reduced pressure approximation displays a greater model discrepancy with higher predicted pressure variations in the reduced model (Figure 5B).

Now, we turn to compare the full and reduced model predictions of physiologically realistic perivascular flow in an image-based PVS surrounding a vascular bifurcation (Model C12). The prescribed pressure difference between inlet and outlets as well as the cardiac wall motion induces pulsatile flow with a net flow component [18] (Figure 6A, Supplementary Video S3). We note that the domain radii vary both angularly and axially, also for the initial domain, and also that the presence of a bifurcation region induces non-Poiseuille/non-Womersley-type velocity profiles. Comparing the full and reduced average pressure and flux at the time of peak velocity (Figures 6F,J), we note that the reduced model captures the qualitative and quantitative flow and pressure characteristics. The bifurcation conditions are satisfied at the bifurcation point b (Figures 6C,G) with a parent branch flux q ̂ ( b ) | Λ P = − 5.52 × 1 0 − 4 μ L/s and daughter branch fluxes q ̂ ( b ) | Λ d 1 = − 7.8 × 1 0 − 5 μ L/s and q ̂ ( b ) | Λ d 2 = − 4.74 × 1 0 − 4 μ L/s. The uneven flux distribution is induced by the smaller average width of one of the daughter vessels. The predicted stress σ ̂ is continuous (data not shown).

The reduced peak cross-section flux (over time) at the inlet is −1.5 × 10⁻³  μL/s, and 1.2 × 10⁻³  μL/s and 7.9, ×, 10⁻⁴  μL/s at the outlets (Figure 6B).

Comparing the relative difference in peak flux at the inlet and outlets, we note that the discrepancy is largest at larger daughter outlet (s _out ) with a relative difference e p p e a k ( s o u t ) = 12.7 % . Comparing the full and reduced peak pressures at the upper (s _u ), middle (s _m ) and lower (s _l ) cross-sections, we find relative differences e p p e a k ( s u ) = 1.1 % , e p p e a k ( s m ) = 4.1 % , and e p p e a k ( s l ) = 16.4 % . The analogous numbers for the fluxes are e q p e a k ( s u ) = 1.0 % , e q p e a k ( s m ) = 33.5 % , and e q p e a k ( s l ) = 5.1 % . Thus, the relative differences in peak flux are larger near the bifurcation region (Figures 6F,J).

The net flow is a key quantity of interest for the physiological relevance of perivascular flow and transport. The net flow per cycle in the full model is 3.5 × 10⁻⁵, and 2.9 × 10⁻⁵ μL for the reduced model, corresponding to a relative difference of 17%.

Accurate direct three-dimensional simulations of pulsatile perivascular fluid flow in large, deforming vascular networks involve a significant computational cost. The expense is dominated by solving large linear systems of equations at each time step. For instance, even the moderate-resolution single-bifurcation model considered here (model C12) includes more than 17 000 vertices, 88 000 mesh cells and 400,000 degrees of freedom. For a small-scale idealized model such as axisymmetric Model A2, the reduced model uses 2.1% of the number of degrees of freedom but approximately the same amount of memory and longer runtime (0.16 vs. 0.35 s per time step, Table 3). However, the one-dimensional models reduce computational cost substantially for the image-based geometries (Table 3). For the image-based perivascular segment (Model B2), the reduced model uses 0.4% of the number of degrees of freedom, 2.0% of the runtime, and 2.1% of the memory of the full model. For the image-based bifurcating PVS (Model C12), the reduced model uses 0.18% of the number of degrees of freedom, 0.6% of the runtime and 2.0% of the memory of the full model. Overall, the reduced model reduces the computational expense, both in terms of computational time and memory, by several orders of magnitude for image-based PVS segments.