A computational framework for quantifying blood flow dynamics across myogenically-active cerebral arterial networks

TL;DR

This paper introduces a new computational method to study how blood flow is regulated in rat cerebral arteries under changing pressure.

Contribution

A novel framework integrating vascular wall mechanics and 1D blood flow dynamics to model cerebral autoregulation in myogenically-active arteries.

Findings

The framework accurately models vascular wall responses to pressure changes at single vessel and network levels.

Weak coupling reduced computational cost while maintaining accuracy for vessel size and boundary conditions.

Pressure surges redistribute flow and vascular constriction across generations of cerebral arteries when myogenic tone is present.

Abstract

Cerebral autoregulation plays a key physiological role by limiting blood flow changes in the face of pressure fluctuations. Although the underlying vascular cellular processes are chemo-mechanically driven, estimating the associated haemodynamic forces in vivo remains extremely difficult and uncertain. In this work, we propose a novel computational methodology for evaluating the blood flow dynamics across networks of myogenically-active cerebral arteries, which can modulate their muscular tone to stabilize flow (and perfusion pressure) as well as to limit vascular intramural stress. The introduced framework integrates a continuum mechanics-based, biologically-motivated model of the rat vascular wall with 1D blood flow dynamics. We investigate the time dependency of the vascular wall response to pressure changes at both single vessel and network levels. The dynamical performance of the vessel wall mechanics model was validated against different pressure protocols and conditions (control and absence of extracellular Ca2+). The robustness of the integrated fluid–structure interaction framework was assessed using different types of inlet signals and numerical settings in an idealized vascular network formed by a middle cerebral artery and its three generations. The proposed in-silico methodology aims to quantify how acute changes in upstream luminal pressure propagate and influence blood flow across a network of rat cerebral arteries. Weak coupling ensured accurate results with a lower computational cost for the vessel size and boundary conditions considered. To complete the analysis, we evaluated the effect of an upstream pressure surge on vascular network haemodynamics in the presence and absence of myogenic tone. This provided a clear quantitative picture of how pressure, flow and vascular constriction are re-distributed across each vessel generation upon inlet pressure changes. This work paves the way for future combined experimental-computational studies aiming to decipher cerebral autoregulation.

Figures

Full text

Introduction

1

Due to their size and extension, small arteries and arterioles are responsible for a significant blood pressure drop across the cerebral circulation (Blanco et al. 2017). Diameter in these resistance vessels is regulated through a combination of local and systemic control mechanisms that operate a muscular apparatus made of smooth muscle cells (SMCs) (Claassen et al. 2021). This enables blood vessels to develop tone across their wall thickness and ultimately to adjust their inner diameter upon different mechanical stimuli. From a hierarchical point of view, the myogenic tone serves as a fundamental, low-level mechanism controlling lumen diameter, as it directly responds to the local pressure level. This mechanism enables the vessel to limit diameter fluctuations and stabilize flow in the face of significant hydrodynamic changes by developing tone within its muscle layer. This means that, in response to an increase in upstream pressure, the vessel will ultimately reduce its diameter in an attempt to maintain a relatively constant flow rate. Haemodynamic forces, such as luminal pressure and shear stress, are sensed in different ways across vascular compartments, having a different impact on vascular contractility (Knot and Nelson 1998; Harraz et al. 2022; Klug et al. 2023). The presence of basal vascular tone is also essential for flow and metabolic control, which operates by relaxing the local vasculature to divert blood flow to the regions where it is needed the most (Longden et al. 2017). From an experimental perspective, isolating these regulatory mechanisms is a cumbersome task, and due to their inter-dependency, the distribution of mechanical forces along the considered vascular segments remains extremely uncertain. Given the key role of mechanical stimuli in cerebral vascular function, developing methods for their systematic quantification is urgently needed.

Computational blood flow dynamics, in conjunction with biologically-motivated vascular wall models, can be used to shed light on different aspects of cerebral autoregulation. Diverse methodologies were proposed to describe the vascular response to acute changes in haemodynamic forces. The seminal work by Carlson and Secomb (2005), Carlson et al. (2008) introduced a high-level framework which integrates together myogenic, flow and metabolic controls without directly including sub-tissue scales, and its wall tension model has been used to investigate the role of myogenic contribution on cerebral blood flow control (Spronck et al. 2012; Zhao et al. 2024). However, in pre-clinical settings, vascular control mechanisms are typically characterized by using a broad repertoire of compounds that selectively activate or inhibit specific cellular components. To elucidate the causative links between drug intervention, luminal pressure and resulting vessel wall deformation, it is essential to develop multi-scale modelling methodologies that can mimic the effect of intracellular biochemical processes on the tissue emergent behaviour. Some authors, including us, have described the pressure-induced vessel contractility in cerebral arteries as a result of intracellular processes driving myosin-actin filament interaction, leading to tone development (Yang et al. 2003a, b; Uhlmann and Balzani 2023; Coccarelli et al. 2024). Yang et al. based their vascular contractility model on a comprehensive description of the pressure-induced $[eqn]$ dynamics in SMC, and assessed the impact of voltage operated $[eqn]$ channel blockage on the static diameter-pressure curve. Conversely, in the work by Uhlmann and Balzani (2023), stretch triggers contraction via a $[eqn]$ -dependent and $[eqn]$ -independent mechanisms, and the model accurately captures the dynamic response of the vascular wall under control or pharmacological inhibition. In our previous work (Coccarelli et al. 2024), we introduced a vascular mechanics model that recapitulates all major pressure-induced intracellular pathways ( $[eqn]$ , ROCK, and PKC) in cerebral arteries. Luminal pressure modulates the activity of distinct intracellular processes, which inform a well-established continuum mechanics-based structural model of active tone generation in SMCs (Murtada et al. 2010; Murtada and Holzapfel 2014; Murtada et al. 2016). These micro-structurally detailed models have been used primarily to study myogenic tone development in isolated cerebral vessels, but there is little quantitative knowledge of how myogenic tone, acting as a feedback mechanism, ultimately redistributes pressure and flow along the vasculature. To our knowledge, only a few studies (Aletti et al. 2016a, b; Daher and Payne 2023) have integrated myogenic contractility into space-dependent blood flow dynamics models of the cerebral circulation. Aletti et al. proposed a 3D fluid dynamic model of retinal circulation, in which the induced vascular structural response is computed using the (mechanical) wall tension model by Yang et al. (2003a), Yang et al. (2003b). On the other hand, Daher and Payne (2023) combined a high-level description of myogenic and shear-dependent regulatory mechanisms with an efficient blood flow model to investigate the impact of autoregulation across a microvascular network. Although they use different methodologies, both studies appear rigorous and report very interesting results. However, since these latter studies do not provide a detailed description of the intracellular processes driving myogenic tone development, their usefulness in modelling the pharmacological modulation of the SMC contractile machinery might have certain limitations.

In this study, we aim to fill this gap by introducing a methodology to integrate a continuum mechanics-based, biologically-motivated vascular wall model within an extensively validated 1D blood flow dynamics framework. This provides a new in-silico tool to predict and recover various haemodynamic and vascular wall quantities under a wide range of conditions, such as variable upstream (or downstream) pressure or the presence of vasoactive agents, which can be replicated in the laboratory.

Methods

2

Blood flow dynamics

2.1

We assume that flow in small cerebral arteries is laminar and axisymmetric, with a Poiseuille velocity profile. The considered 1D fluid domain may range from a single vessel to a complex vessel network. The mass and momentum conservation equations for a fluid flowing in a collapsible vessel can be written in the pressure-flow form as Carson and Van Loon (2017), Coccarelli et al. (2021):

$$[eqn]$$

where $[eqn]$ is the axial direction, $[eqn]$ is the luminal cross-sectional area, $[eqn]$ is the average pressure in the cross section corresponding to pressure acting on the inner wall surface, $[eqn]$ is the volumetric flow rate in the cross section whilst $[eqn]$ and $[eqn]$ are respectively the blood density and dynamic viscosity which, for the sake of simplicity, are assumed constant. It is worth mentioning that the variation of the area with respect to the fluid pressure defines the vessel compliance $[eqn]$ and can be determined from the constitutive law of the vascular wall. In line with Carson and Van Loon (2017), the system (1) is linearized with respect to time as follows

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the current and previous time steps. The fluid domain is subdivided into elements of non-necessarily equal size. Following the work by Carson and Van Loon (2017), Eqs. (2) are integrated in space using the enhanced trapezoidal rule method and discretized in time using a second-order backward difference scheme. After some steps, the system of equations in (2) may be re-written at the element level in the following compact form:

$$[eqn]$$

in which $[eqn]$ represents the elemental level, $[eqn]$ , $[eqn]$ are the stiffness matrices of pressure and flow, $[eqn]$ and $[eqn]$ the vectors containing the current element values of pressure and flow, and $[eqn]$ the vector representing convection and diffusion components evaluated at previous time step. Equation (3) serves for the assembling of the global system matrix, which, in conjunction with the boundary conditions, are used to compute (through ‘spsolve’ from SciPy 1.6.0) the pressure and flow rate at the next time step ( $[eqn]$ ) across the whole fluid domain. Bifurcations, unifications, vessel geometry, and material discontinuities are handled in line with Carson and Van Loon (2017). The luminal cross-sectional area and wall compliance at the next time step ( $[eqn]$ ) are subsequently recovered from the vascular wall constitutive relationship, as reported in the following.

Vascular wall dynamics

2.2

The mechanics across the vascular wall are described by following the bottom-up approach reported in our previous work (Coccarelli et al. 2024). To describe pressure-induced tone generation, this model was developed and calibrated using ex vivo data from rat middle cerebral arteries, with the endothelium and nerve contributions removed (Johnson et al. 2009). The wall of cerebral vessels is a complex structure endowed with the capacity to generate tone upon luminal pressure loading. Although the wall consists of functionally distinct layers, we assume that its entire volume is occupied by SMCs, which are partitioned into $[eqn]$ homogeneous cellular domains ( $[eqn]$ ) along the vessel’s radial direction. The contractile activity of SMCs is described by considering (luminal) pressure-induced biochemical signalling, which modulates myosin-actin interactions alongside cytoskeleton remodelling. The intracellular information serves as input for the structural mechanics model, which determines the emergent tissue response to pressure.

Intracellular chemo-mechanics

2.2.1

The pathways activated by pressure in SMCs are represented through the following (normalized) quantities: $[eqn]$ concentration $[eqn]$ , ROCK activity level $[eqn]$ , HSP27 phosphorylation level ( $[eqn]$ ), MLCP phosphorylation level ( $[eqn]$ ), Cofilin phosphorylation level $[eqn]$ phosphorylation level $[eqn]$ , and G-actin content ( $[eqn]$ ). The time evolution of these variables is evaluated through a logic-based signalling graph, which is translated in the following system of equations:

$$[eqn]$$

where $[eqn]$ with $[eqn]$ is a logistic function connecting two signalling variables whilst $[eqn]$ with $[eqn]$ represents the time constant of the associated intracellular process. Equations (4) allow us to evaluate how changes in luminal pressure (over time) influence cross-bridges (XBs) formation (between actin and myosin) and cytoskeleton remodelling, which are represented by $[eqn]$ and $[eqn]$ , respectively. These molecular factors and SMC stretch level drive the relative sliding between actin and myosin filaments, which, together with cytoskeleton stiffness (represented by the F-actin content), enable tone development. Given a luminal pressure level $[eqn]$ , Eqs. (4) are discretized in time with the two-step Adams-Bashforth method and solved through a nonlinear equations solver (‘root’ from SciPy 1.6.0, with default settings). For the sake of simplicity, all the SMCs are assumed to be aligned along the circumferential direction, and only the corresponding stretch component $[eqn]$ plays a significant role in tone development. The cell contractile fibres (CFs) can be represented as a series of interconnected contractile units (CUs) that are anchored at the cell membrane through actin cortex passive elements (see (Murtada et al. 2016) for more details). The dynamics governing the (normalized) relative filament sliding ( $[eqn]$ ) within each CU is expressed via

$$[eqn]$$

where $[eqn]$ is the time constant associated with actin-myosin filaments sliding dynamics (and force generation). $[eqn]$ is the average driving force generated from the XBs cycling

$$[eqn]$$

where $[eqn]$ is the average length of myosin filaments, $[eqn]$ is the maximum phosphorylation rate, $[eqn]$ is the average distance between myosin monomers heads, $[eqn]$ represents the XB elastic stiffness and $[eqn]$ is the average displacement associated to power-stroke, whilst $[eqn]$ describes the filament overlap, which depends on the relative filament sliding $[eqn]$ via

$$[eqn]$$

The reaction force due to the resistance from the number of contractile units ( $[eqn]$ ) in series with a F-actin element at each extremity is given as the product between the total (XBs and passive elements) elongation and the resulting stiffness of the contractile fibres:

$$[eqn]$$

where $[eqn]$ is the stiffness associated with a number of CUs $[eqn]$ , which is directly related to the level of LC_20_ phosphorylation via

$$[eqn]$$

whilst $[eqn]$ is the actin cortex stiffness, which can be evaluated as a function of the F-actin content level

$$[eqn]$$

where $[eqn]$ is the maximum stiffness under loading conditions, whilst $[eqn]$ and $[eqn]$ are the coefficients of the associated activation function.

Tissue mechanics

2.2.2

The vessel wall consists of multiple concentric SMC layers, and here, in line with Coccarelli et al. (2021), Coccarelli et al. (2024), it is modelled as an axisymmetric homogeneous hyperelastic thick-walled tube. The stretch $[eqn]$ experienced by an SMCs domain coincides with the circumferential stretch of that portion of tissue. Each cellular domain is governed by a distinct intracellular chemomechanics, which depends on the average $[eqn]$ within the domain. For the sake of simplicity, the biochemical signalling component depends only on luminal pressure and therefore does not vary between cellular domains. However, the model’s structure allows for the future incorporation of space-dependent stress levels as input. Finite strain theory and incompressibility assumption are adopted for describing the tissue kinematics

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the radial and axial stretches, respectively, whilst $[eqn]$ is the reference radius, $[eqn]$ the deformed radius and $[eqn]$ a parameter accounting for the residual strain. The usual mapping between reference ( $[eqn]$ ) and deformed ( $[eqn]$ ) configuration is used:

$$[eqn]$$

where $[eqn]$ is the deformed luminal radius, whilst $[eqn]$ and $[eqn]$ are the luminal radius and thickness in the reference configuration, respectively. The luminal circumferential stretch is indicated with $[eqn]$ while the deformed luminal and outer diameters are indicated by $[eqn]$ and $[eqn]$ , respectively. The generated tone from CFs can be quantified through the average first Piola–Kirchhoff stress $[eqn]$ of the vascular tissue, $[eqn]$ , where $[eqn]$ is the surface density of the CFs. The momentum conservation principle along the radial direction allows us to link the luminal cross-sectional area ( $[eqn]$ ) to luminal pressure $[eqn]$ (Coccarelli et al. 2021):

$$[eqn]$$

where $[eqn]$ is the external pressure acting on the outer surface of the vessel, $[eqn]$ is the material strain-energy function, which is decomposed into active and passive components. The former can be evaluated as

$$[eqn]$$

whilst the passive behaviour is described in line with Holzapfel et al. (2000)

$$[eqn]$$

where $[eqn]$ , $[eqn]$ and $[eqn]$ are the media constitutive parameters, $[eqn]$ , $[eqn]$ with $[eqn]$ being the orientation angle of a collagen fibres family, which are oriented along the circumferential direction of the vessel. While Eq. (13) concerns the momentum balance across the whole wall thickness, Eq. (5) describes the time dependency of quantities that depend on the radial coordinate $[eqn]$ . The integral in Eq. (13) depends on both $[eqn]$ and $[eqn]$ since $[eqn]$ and is computed via Simpson’s rule. In this study, the adopted number of integration points (five) coincides with the number of cellular domain $[eqn]$ and therefore $[eqn]$ and $[eqn]$ are computed at the middle of each cellular domain and represented by $[eqn]$ , $[eqn]$ with $[eqn]$ . To evaluate the mechanics across the vascular wall we rewrite Eqs. (5,13) as residuals

$$[eqn]$$

The set of Eqs. (16) is discretized in time with the two-step Adams–Bashforth method and solved by using a solver for non-linear equations (‘root’ from SciPy 1.6.0, tol=1e-8) providing the solution ( $[eqn]$ ) at the next time step. Once the area $[eqn]$ is known, the compliance $[eqn]$ can be evaluated using a first-order centred finite difference scheme as done in our previous work (Coccarelli et al. 2021). Alternatively, the vessel wall mechanics can be evaluated by assuming that the active stress component depends only on the average trans-mural $[eqn]$ , which is obtained by considering the average circumferential stretch over the vessel thickness. This strategy, already employed in Coccarelli et al. (2024), is more computationally efficient than the original without ‘averaged active stress’ as it requires simultaneously solving only two equations rather than the set of Eqs. (16).

Multi-physics coupling strategy

2.3

The fluid and solid wall equations are strongly coupled via the fixed-point iteration method (Fig. 1). The guess values for the cross-sectional area $[eqn]$ across the whole vascular network are initialized by considering the luminal areas at the previous time step. The numerical solution of 1D fluid flow equations (reported in 2.1) provides the pressure and volumetric flow fields ( $[eqn]$ ) throughout the fluid domain. Each fluid node is associated with a vascular tissue ring, whose mechanics does not depend on the neighbouring nodes along the vessel’s axial direction. The vascular chemo-mechanical model (reported in Sect. 2.2) is used to evaluate sequentially the intracellular signalling variables $[eqn]$ and $[eqn]$ and the resulting luminal area $[eqn]$ , alongside the relative filament sliding $[eqn]$ . In the current methodology Eqs. (4) only depend on the current luminal fluid pressure $[eqn]$ and, therefore, do not need to be solved simultaneously with Eqs. (16). The new value for the cross-sectional area $[eqn]$ is used to update the wall compliance $[eqn]$ . If the Root Mean Square Relative Error (RMSRE) between $[eqn]$ and $[eqn]$ across the vascular network is below the prescribed tolerance $[eqn]$ , the variables $[eqn]$ , $[eqn]$ , $[eqn]$ are stored as the solution of the current time step; otherwise the fluid-structure interaction sequence is again re-iterated with $[eqn]$ as the new guess values.

The tolerance for the fluid–structure interaction coupling ( $[eqn]$ ) is set to 0.000001 throughout the study, unless specified otherwise. In the following, we refer to ‘strong coupling’ when the solution is obtained by applying the fixed-point iteration method, while to ‘weak coupling’ when fluid and solid equations are solved sequentially, but there is no iterative procedure for updating the luminal area. In this study, the entire model was implemented in Python 3.8, and all the simulation results were obtained using a desktop workstation with an Intel(R) Core(TM) i7–9700K CPU @ 3.60GHz.

Time constants and variables initialization

2.4

The behaviour of the vascular contractility model was investigated at steady state in our previous work (Coccarelli et al. 2024) under different conditions (control and selective pharmacological modulation). By considering a cannulated vessel (no flow, same pressure at the extremities) immersed in a bath of physiological saline solution (PSS), here we aim to characterize its time-dependent response to upstream pressure variations. The dynamics of the SMC signalling are governed by the time constants ( $[eqn]$ ), each of them reflecting the speed of variation of an intracellular process. Reference experimental studies (Johnson et al. 2009; Moreno-Dominguez et al. 2014) reported diameter-vs-time recordings of different pressurized arteries following the addition of ROCK and PKC inhibitors to the PSS bath. We observed that the duration of the transient varied significantly across the reported dataset, and the number of traces may be too small to fully characterize the dynamics of the processes. Furthermore, the diffusion kinetics of the drug across the bath represents another source of uncertainty. In this study, we assume that the total time required for each individual pressure-induced pathway ( $[eqn]$ , ROCK, PKC) to convert the mechanical stimulus (luminal pressure) into a signal for the contractile apparatus (actin-myosin filaments and cytoskeleton) is the same. This implies that the transduction of the applied load information into filament sliding and cytoskeleton remodelling is governed by a single time constant, $[eqn]$ , while $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ are considered infinitesimal. An early study (Kontos et al. 1978) suggests that the diameter response to pressure variation in cerebral arteries develops within a fraction of a minute. Therefore, for this time constant, we initially considered different values: $[eqn]$ , and 60 s. The dynamics governing the actin-myosin filaments sliding within the passive surrounding matrix is characterized by the time constant $[eqn]$ , which was estimated in the study by Murtada et al. (2016) for mouse descending aorta $[eqn]$ s). Due to the different vessel size and functionality between the latter and cerebral arteries, we also considered alternative values around the one proposed by the previous study ( $[eqn]$ ). In all numerical experiments conducted in this work, the intracellular variables are initialized to 0, except for $[eqn]$ which is set to −0.02, while the initial $[eqn]$ and $[eqn]$ are set equal to the load-free values.

Results

3

Pressure-induced arterial wall dynamics: comparison versus experiments

3.1

Alternated pressure protocol

3.1.1

The dynamical behaviour of the vascular structure is assessed by simulating its response to alternated luminal pressure levels (between 10 and 60 mmHg) under control (PSS) and pharmacological modulation of the $[eqn]$ pathway (via the addition of 30 μM Diltiazem to the PSS). The effect of $[eqn]$ on all the modelled intracellular quantities upon pressure-activation is reported in the Appendix. Pressure variations within the 10–60 mmHg range induce a remarkable change in wall $[eqn]$ concentration, while variations in ROCK activity and consequent MLCP phosphorylation are limited. This is in line with the hypothesis for which some $[eqn]$ sensitization mechanisms are mainly activated at medium-high pressure levels (> 60 mmHg) (Osol et al. 2002). On the other hand, it is difficult to assess and verify the simulated dynamics of the PKC pathway due to its unresolved bidirectional dependency with $[eqn]$ activity in cerebral SMCs (Osol et al. 1991; Gokina et al. 1999; Earley et al. 2007; El-Yazbi et al. 2015). Model predictions are compared against the experimental recordings in rat distal posterior cerebral arteries by Knot and Nelson (1998). In this numerical experiment, we only use the vascular mechanics model (described in Sect. 2.2) as the wall deformation and relative filament sliding are calculated by prescribing the luminal pressure as a function of time (see the Appendix for more details on the signal) and no ‘averaged active stress’ is adopted. From Knot and Nelson (1998), we assume the load-free ( $[eqn]$ ) outer diameter for the vessels in the control and inhibited case equal to 135 and 150 μm, respectively. All the other parameters of the wall mechanics model remain the same as in our previous study (Coccarelli et al. 2024) (here and in the following sections), unless specified otherwise. For this simulation, the time step is set to 0.5 s. Diltiazem at a concentration of 30 μM is expected to prevent the opening of L-type calcium channels upon pressure-induced membrane depolarization, restricting the $[eqn]$ influx into the intracellular space, which causes a reduction in tone development. Here, we assume that, upon drug effect, pressure has a limited effect on the intracellular $[eqn]$ concentration, with only a marginal increase due to stretch-operated channels (see the Appendix for more details on the adopted pressure - $[eqn]$ relationship). By considering different combinations of $[eqn]$ and $[eqn]$ (as described above), we estimate the associated RMSRE between the simulated and experimental outer diameter traces under control condition. The combination of parameters ( $[eqn]$ ) yields the smallest discrepancy between simulated and experimentally-recorded outer diameter (Table 1).

The RMSRE seems to be more sensitive to changes in the intracellular signalling time constant $[eqn]$ , rather than to the time constant associated with the actin-filaments sliding $[eqn]$ . Figure 2 shows how the intracellular $[eqn]$ concentration and the outer diameter change in cerebral arteries when alternated luminal pressure levels are applied for the control and inhibition (30 μM Diltiazem) cases.

As expected, the timescale of the pressure-induced intracellular processes ( $[eqn]$ ) plays a profound effect on the shape of the cytosolic $[eqn]$ concentration as well as on the outer diameter responses. In the control case, a smaller $[eqn]$ corresponds to a higher $[eqn]$ increase/decrease rate, while larger time constants are associated with a significantly slower variation in $[eqn]$ . Although $[eqn]$ closely matches the $[eqn]$ drop after pressure lowering, smaller time constants better represent the pressure-induced $[eqn]$ elevation. The disparity in optimal $[eqn]$ values between the pressure-induced $[eqn]$ increase and decrease may be due to the different $[eqn]$ dynamics associated with each phase. The effect of different combinations $[eqn]$ on the intracellular $[eqn]$ concentration and outer diameter responses in the control case is reported in the Appendix. Among the considered values, $[eqn]$ is the best choice for capturing the experimental $[eqn]$ and diameter traces over time in the control case, and this closely aligns with previous research (Daher and Payne 2023).

Multiple steps pressure protocol

3.1.2

The dynamical response of the vascular wall to pressure variations is expected to depend also on the geometrical features of the vessel itself, such as the wall thickness. To explore this aspect and further validate our vessel wall dynamics model, we consider another experimental protocol in which the inflating pressure is increased from 10 to 120 mmHg through multiple sequential steps (see Fig. 3).

We simulate the changes in diameter over time for three levels of the load-free thickness-medium radius ratio $[eqn]$ in the control and 0 extracellular $[eqn]$ cases, and compare them against representative experimental time recordings of diameter from rat middle cerebral arteries (Johnson et al. 2009). For this simulation, the load-free diameter for the control and $[eqn]$ removal cases are 190 and 230 μm, respectively, while the time step is set to 0.01 s. The experimental traces obtained under control conditions show two distinct qualitative behaviours while the $[eqn]$ removal cases share a much more similar trend. Our proposed vascular wall model predicts diameter variations within the same range as the measurements, for both control and 0 extracellular $[eqn]$ conditions. In both cases, a thicker vessel generates more tone and limits diameter variation despite increasing pressure. We speculate that the difference in behaviour between the experimental control traces may be explained, at least in part, by an effect of vessel size.

Pressure-induced dynamics across an arterial network

3.2

The effect of myogenic tone on blood flow regulation can be well appreciated by evaluating haemodynamic quantities across a cerebral arterial network. Since the focus of this work is on the definition of a suitable methodology for modelling blood flow within self-regulated vessels, we consider an idealized symmetrical network branching from a rat’s middle cerebral artery. Here, we introduce a vascular network (morphology reported in Fig. 4) to assess i) the accuracy and efficiency of the different solution procedures and ii) the effect of upstream pressure changes on the system’s blood flow dynamics. The generation $[eqn]$ represents the last artery that precedes the arteriolar vasculature in the parenchymal space. Parenchymal arterioles also develop myogenic tone, and they are expected to play an important role in the stabilization of flow and perfusion pressure in the face of upstream pressure changes (Iadecola 2017). However, their structure and intracellular signalling factors may differ from cerebral arteries (Cipolla et al. 2014; Li and Brayden 2017). For the sake of simplicity, we do not explicitly represent these vascular beds (and their downstream vessels) and the network is truncated after the $[eqn]$ vessels. The downstream circulation is represented through an (outlet) pressure boundary condition $[eqn]$ , and backward pressure wave reflections are prevented by including a characteristic impedance $[eqn]$ between the terminal node of each $[eqn]$ vessel and the associated outlet node (where $[eqn]$ is imposed). To reflect in vivo conditions, $[eqn]$ and $[eqn]$ are set to 50 mmHg and 10 mmHg, respectively.

Blood viscosity and density are set to 0.05 poise and 1.04 g · cm^−3^, respectively. In each vessel, the axial spatial domain is discretized with two elements. Here and in the following, the pressure is initialized with the first inlet value whilst the initial flow rate is set to 0 ml/s.

Comparison between solution procedures

3.2.1

Here, we evaluate the impact of numerical procedure settings such as the time step $[eqn]$ , coupling type and active stress averaging on the accuracy and efficiency of the solution across the network. To appropriately test the model robustness, we prescribe a periodic pressure signal with a mean that varies over time at the inlet. We consider the simulation results obtained with $[eqn]$ , strong coupling and without averaged active stress as the ground truth solution (case 1).

Dramatic changes in upstream pressure are accompanied by substantial variations in flow rate, luminal diameter and actin-myosin filament sliding across the vessel network (Fig. 5). As expected, the fluctuation amplitude in all the recorded variables is significantly more mitigated in the higher-generation vessels than in the upstream larger arteries. Alternatively, more computationally efficient simulation settings than case 1 are explored (Table 2). The numerical accuracy of cases 2–5 was assessed (against case 1) by evaluating how RMSRE for flow ( $[eqn]$ ) and luminal area ( $[eqn]$ ) distribute across the vascular network (Table 3).

Simulation results indicate that a time step equal to 2.5 × 10^−4^ s provides an optimal compromise between numerical accuracy and computational speed. Weak coupling provides another significant reduction in Wall Clock Time (WCT) without remarkably affecting the solution precision (all $[eqn]$ RMSREs in case 4 are actually slightly improved vs case 3). This information may be extremely important when the model is used to describe vast vascular networks and/or simulate more comprehensive cellular dynamics. Adopting an ‘averaged active stress’ enables a further significant reduction in computational time and maintains the error below a reasonable threshold (Mean $[eqn]$ RMSRE < 1.9 %, Mean $[eqn]$ RMSRE < 0.16 %). The discrepancy between the numerical solutions obtained with the most and least computationally demanding strategies (cases 1 and 5) is shown in Fig. 6.

The time evolution of the relative error for flow rate and luminal area of cases 4 and 5 (with respect to case 1) across the vascular network indicate that in both cases the accuracy of the solution does not significantly deteriorate over the considered time interval (see Appendix for the error comparison between cases). The numerical strategy associated with case 5 is adopted for the simulations in the next section.

Myogenic response to upstream pressure increase

3.2.2

When the upstream pressure increases, the myogenic tone enables small arteries (and arterioles) to adjust their luminal diameters to stabilize blood flow and limit perfusion pressure variations. Figure 7 shows the myogenic response across all the vessel generations to an extreme upstream pressure change from 50 to 120 mmHg.

The pressure surge initially causes an increase in diameter and flow, which are then gradually reduced as the tone develops. The propagated pressure change from upstream is enormously mitigated at higher-generation vessels. Alongside the curves representing myogenically-active arteries (control), we also reported the predictions for vessels with impaired contractile capacities due to extracellular $[eqn]$ removal (0 $[eqn]$ ). The comparison between these vessel conditions highlights the importance of myogenic tone in counteracting acute hydrodynamic changes. Under control conditions, all vessels constrict to redistribute the new pressure load and minimize flow variations, while in the presence of tone (partial) inhibition, the vessels dilate upon pressure increase. The variation of $[eqn]$ in time shows how functioning contractile units respond to a pressure surge. We evaluated the impact of the final (steady state) inlet pressure level on the total flow rate through the vascular system (recorded at $[eqn]$ ). The steady-state results obtained in the control and 0 $[eqn]$ cases define two distinct relationships between upstream pressure and flow (Fig. 8). Upon inlet pressure increase, control conditions are associated with moderate flow rate increments, whilst the 0 $[eqn]$ case exhibits almost an exponential dependency between flow rate and inlet pressure.

The luminal diameter ratio (control over 0 $[eqn]$ ) is used to define the level of vessel constriction across the network for various final upstream pressures. For the boundary conditions considered, large myogenically-active vessels are subjected to higher luminal pressure and provide more resistance to flow than the smaller arteries. Although the proposed numerical experiments allow us to quantify the impact of myogenic tone on blood flow stabilization across a network of (rat) cerebral arteries, the associated results may vary depending on the network morphology and the role of the downstream circulation.

Discussion

4

Thanks to myogenic tone, cerebral arteries and arterioles can locally adjust their diameter to maintain nearly constant blood flow and perfusion pressure despite acute changes in upstream pressure. Therefore, this local blood flow control helps prevent eventual local tissue oxygen and nutrient starvation. Malfunctioning of this regulatory mechanism can contribute to the development and progression of different life-threatening conditions, and indeed, its restoration has been proposed as a potential therapeutic target (Hill and Al 2009; Palomares and Cipolla 2014; El-Yazbi and Abd-Elrahman 2017; Lidington et al. 2018). Despite its physiological importance and implication in diverse pathological conditions, the effect of myogenic tone on blood flow dynamics in cerebral vessels has been rarely analysed and quantified through in-silico approaches. The few currently available models, although computationally efficient, do not directly relate the intracellular processes underlying myogenic tone to blood flow changes across vascular networks. Therefore, there is an urgent need for more physiologically-representative modelling methodologies capable of dealing with different types of data from the laboratory and clinic (Payne 2024).

Through this work, we propose a comprehensive computational methodology which integrates, in a robust manner, the vascular SMC contractile machinery into blood flow dynamics in cerebral arteries. This framework aims to complement, confirm and conceive future investigations on blood flow regulation due to myogenic tone across the cerebral circulation. The time-dependent behaviour of the model was investigated at both single vessel and network levels. The comparison with the time recordings of the vessel wall $[eqn]$ concentration and diameter from the alternated pressure protocol allowed us to identify the time constants governing the dynamical response. The results obtained from the simulated drug intervention (30 μM Diltiazem, alternated luminal pressure protocol) and the diameter predictions under increasing pressure via multiple sequential steps further validated the proposed methodology. Regarding the second pressure protocol, the experimental traces reported in Johnson et al. (2009) exhibit noticeable variability in behaviour with increasing pressure. This variability may be due to several factors including axial stretch (Bell et al. 2017) or simply differences in vessel size. Our model demonstrates how the latter can dramatically influence the wall’s response to increasing pressure. Altogether, these results demonstrated that the developed model can accurately capture the pressure-induced wall dynamics under different conditions (RMSRE < 14 % in the control diameter under alternated luminal pressure using the identified time constants).

The main aim of this work was to assess the proposed fluid–structure interaction methodology in terms of numerical convergence and performance. To do this, we defined an idealized network made of a rat’s middle cerebral arteries and its three (symmetric) generations. The simple morphology of the vascular network hindered the comparison against experimental flow measurements but allowed to clearly quantify the haemodynamics for each vessel category. We prescribed a highly variable signal for the pressure at the inlet whilst reflection-free conditions were imposed at the outlets of the vascular network. By considering this newly defined toy-problem, we analysed the impact of different numerical settings on the computational speed and accuracy of the solution. Weak coupling and considering an average active stress allows for a substantial reduction in computational time, without significantly sacrificing the accuracy of the solution. Imposing the reflection-free outflow conditions ignores backward components from the downstream circulation but allows for an exclusive assessment of the $[eqn]$ vessels’ contribution to the network flow dynamics. Future research on the impact of outflow conditions on vascular network dynamics is warranted. Overall, the numerical strategy of case 5 represents a good compromise between accuracy and computational costs for problems with a time span of a few minutes. For problems with longer time duration and where high accuracy is necessary, we recommend weak coupling without average active stress.

The model was then used to assess how a change in pressure at the inlet of a middle cerebral artery is mitigated (via myogenic mechanism) across an arterial network composed of its three generations. The computed results highlight the importance of myogenic tone in limiting blood flow variation upon significant pressure changes. The proposed framework is based on a biologically-motivated model for the vascular wall (Coccarelli et al. 2024) whose parameters were identified by considering both intracellular and tissue recordings under different conditions (control, no extracellular $[eqn]$ , with vasoactive agents). Since the role of the endothelium on myogenic tone is not yet fully elucidated in cerebral resistance arteries (McCarron et al. 1989; Harder 1987; Wallis et al. 1996), we have not explored it here but plan to do so in a future study. To simplify the analysis of the proposed model, we have assumed that the $[eqn]$ , ROCK and PKC pathways are governed by the same time constant. This simplification may need to be reconsidered once quantitative experimental data on the timing of distinct intracellular processes become available. With these settings, the current framework predicts that any increase in upstream pressure will lead to a moderate rise in blood flow through the arterial network (compared to the case without extracellular $[eqn]$ ). This might seem to slightly contrast with the ‘theoretical’ autoregulation curve ‘flow rate - upstream pressure’ (see (Aletti et al. 2016a) for instance), for which flow remains constant across a mid-pressure interval. However, previous studies (Bryan et al. 2001; Toth et al. 2011) have demonstrated that luminal flow also contributes to tone development in cerebral arteries, in a manner distinct from that of large extra-cranial arteries (Carter et al. 2016; Smith et al. 2019). In pressurized pial vessels, an increase in flow causes further vasoconstriction, and this seems to synergistically work together with the myogenic mechanism to preserve blood volume within the intracranial space (Koller and Toth 2012). The current modelling framework does not account yet for the flow-induced tone regulation, and this may explain why the predicted flow rate is not maintained constant across the upstream pressure range 60–100 mmHg but presents a linear dependency that is consistent with the data reported in Toth et al. (2011). We are currently investigating this regulatory component and its integration into our computational framework. Different modelling approaches are available today to quantify and characterize myogenic regulation in the cerebrovasculature. Modelling SMC contraction in 3D tissue samples (Murtada and Holzapfel 2014; Coccarelli et al. 2018; Uhlmann and Balzani 2023) allows to explore the stress field across the wall with high resolution but its application to vessel networks is hindered by its exorbitant computational cost. To enhance the model’s capability and reliability in simulating pharmacological interventions, more detailed representations of the underlying SMC processes (Yang et al. 2003a, b; Kenny et al. 2018; Hernandez-Hernandez et al. 2024) could be integrated into the framework in the future. However, our reduced-order model still requires several hours of computational time to simulate the myogenic response for a few minutes in a network of small rat arteries. Due to the small geometrical size of the considered vascular network, a severe time step was employed to maintain a low Courant–Friedrichs–Lewy (CFL) number. Various strategies, including dual time-stepping and machine learning, can be explored in the future to enhance the framework’s scalability for larger vascular networks or to incorporate new cellular regulatory mechanisms (Kenny et al. 2018; Sten et al. 2023). Combining active vascular wall models with compartmental/0-D blood flow models (Ursino and Giannessi 2010; Spronck et al. 2012; Tong et al. 2021; Sten et al. 2023) is an attractive alternative due to its efficiency and ability to provide a holistic view of integrated autoregulatory mechanisms. In certain situations, using fewer parameters may also facilitate the integration of experimental recordings into the model. On the other hand, their application is limited when it comes to quantifying the spatial distribution of haemodynamic forces and vessel diameters along a vascular network, as well as simulating pharmacological interventions. Overall, we believe that the proposed modelling framework offers a balanced trade-off between detail accuracy, computational efficiency, and scalability to relevant pre-clinical scenarios.

To conclude, the presented framework represents an essential tool for investigating and quantifying blood flow regulation mechanisms in cerebral circulation. If provided with flow information at the inlet and outlet, the model can recover the mechanical stimuli acting on the SMCs along the vascular network, as well as their contribution to blood flow regulation. While this model can already be used to mimic some ex vivo scenarios, we plan to extend its predictive capacity by incorporating new modelling components (such as flow-induced tone and metabolic function) and validating them against new in vivo data. This framework serves as a tool to quantify the role of myogenic tone in cerebral blood flow and ultimately aims to improve the accuracy of computer predictions in clinical scenarios (Sun et al. 2024).