A Multicellular Vascular Model of the Renal Myogenic Response

Abstract

The myogenic response is a key autoregulatory mechanism in the mammalian kidney. Triggered by blood pressure perturbations, it is well established that the myogenic response is initiated in the renal afferent arteriole and mediated by alterations in muscle tone and vascular diameter that counterbalance hemodynamic perturbations. The entire process involves several subcellular, cellular, and vascular mechanisms whose interactions remain poorly understood. Here, we model and investigate the myogenic response of a multicellular segment of an afferent arteriole. Extending existing work, we focus on providing an accurate—but still computationally tractable—representation of the coupling among the involved levels. For individual muscle cells, we include detailed Ca${}^{2+}$ signaling, transmembrane transport of ions, kinetics of myosin light chain phosphorylation, and contraction mechanics. Intercellular interactions are mediated by gap junctions between muscle or endothelial cells. Additional interactions are mediated by hemodynamics. Simulations of time-independent pressure changes reveal regular vasoresponses throughout the model segment and stabilization of a physiological range of blood pressures (80–180 mmHg) in agreement with other modeling and experimental studies that assess steady autoregulation. Simulations of time-dependent perturbations reveal irregular vasoresponses and complex dynamics that may contribute to the complexity of dynamic autoregulation observed in vivo. The ability of the developed model to represent the myogenic response in a multiscale and realistic fashion, under feasible computational load, suggests that it can be incorporated as a key component into larger models of integrated renal hemodynamic regulation.

1. Introduction

The kidneys are organs in the body of mammals that extract blood waste and regulate the balance of water, electrolytes, acid-base species, etc. [1]. These tasks are crucial for the whole mammalian body and are accomplished by complex processes involving, among others, filtration and reabsorption mechanisms that take place exclusively within the kidney’s nephrons [1]. For these mechanisms to operate properly, the rate of blood volume that is filtered into each of the nephrons must be maintained within a narrow range [1]. Indeed, if the filtration rate is high, nephrons may not have sufficient reabsorbtive capacity and necessary blood substances may be lost in urine. On the other hand, if the filtration rate is low, unnecessary substances or toxic waste may be retained and leak back into the general circulation.

To maintain a stable filtration rate under physiologic conditions, blood supply to the nephrons is regulated individually by a number of autoregulatory mechanisms [2,3,4,5]. All such mechanisms that have been identified to-date share a common effector, namely the afferent arteriole, which is the vessel supplying blood to each nephron [5].

One of the strongest mechanisms acting on the afferent arteriole is the myogenic response [6,7]. This mechanism induces vasoconstriction of the arteriole when blood pressure is elevated and vasodilation when blood pressure is reduced. This way the myogenic response enables the kidney to actively impede or accelerate blood flow and therefore to adjust the filtration rate in response to perturbations of blood pressure that otherwise have huge destabilizing effects [7,8,9,10].

At the organ or vascular levels, the myogenic response in the kidney has been well characterized experimentally in vitro as well as in vivo [11,12,13,14,15,16,17,18]. However, theoretically a deeper understanding of the involved processes is lacking which limits our ability to investigate in silico realistic scenarios under health or disease conditions. For example, although several models accurately represent the myogenic response at the vascular or supravascular levels [2,7,11,19,20,21,22,23,24,25,26,27,28], and few even at the cellular level [6,29,30], none of the existing models attempt to model the myogenic response across multiple scales.

Indeed, multiscale modeling that spans processes ranging from the subcellular up to the organ level is required to elucidate clinical and experimental findings such as the effects of integrins [31,32], epithelial sodium channels (ENaCs) [33,34], angiotensin II (AngII) [35,36,37,38], nitric oxide (NO) [39,40], and an array of other vasomodulators [41,42,43,44], or the pathophysiology associated with certain chronic conditions such as hypertension [45,46], diabetes mellitus [47,48], or chronic kidney disease [49,50]. Most of these previous studies have used simple compartmental models (Windkessels), ordinary or partial differential equations, and more recently, probabilistic methods for reproducing renal vascular networks [28].

Our previous modeling study [29] considered only a single afferent arteriole smooth muscle cell and investigated the myogenic response initiated by pressure changes. In [7], we considered a model exhibiting the myogenic response and consisting of an afferent arteriole segment, glomerular filtration, and a short loop of Henle. The smooth muscle cell in this study is less realistic than the model in [29] since it includes a small number of components, such as membrane potential and calcium concentration. In [30], we extended the model in [29] to a segment of multiple smooth muscle cells connected through gap junctions; one limitation of this model is that the numerical implementation (fractional splitting) did not allow us to consider a realistic cell size.

A goal of this study is to gain a better understanding of the interactions developed between subcellular, cellular, and vascular processes involved in the initiation and development of the myogenic response. Another goal is to formulate a realistic model entailing a feasible computational cost that may be further extended to larger scales. To that end, we expand on our existing, highly detailed mathematical model of Ca${}^{2+}$ signaling within an individual afferent arteriole smooth muscle cell [29] that represents the transmembrane transport of major cytosolic ions, intracellular Ca${}^{2+}$ dynamics, the kinetics of myosin light chain phosphorylation, and the mechanical behavior of the entire cell and local vascular wall. We extend this model to a multi-cell vascular model of the afferent arteriole and couple it with a representation of hemodynamics, i.e., blood flow and pressure, and intercellular conduction that also incorporates a layer of endothelial cells. Our model consists of a mixture of ordinary and partial differential equations and we apply numerical methods for differential-algebraic systems. We use the resulting model to study myogenically induced responses of the entire afferent arteriole to characteristic time-independent and also time-dependent pressure perturbations.

2. Mathematical Model

In this section we describe the model equations, the numerical methods developed for the solution of the model equations, and present values used for the model parameters.

2.1. Model Description

To represent the processes mediating the initiation and progress of the myogenic response in the renal afferent arteriole, we model a segment of the mammalian renal vasculature consisting of multiple cells. We incorporate a simplified geometry (Figure 1) and dimensions consistent with the afferent arterioles associated with the long nephrons of the rat kidney [51,52].

To obtain reliable predictions of the driving signal of the myogenic response, i.e., blood pressure profile along the model vessel [7], we represent hemodynamics. Vascular diameter across the model vessel greatly impacts hemodynamics [7], and we combine the myogenic mechanisms of a large number of smooth muscle cells that are the main determinants of vascular tone and wall mechanics [3]. For the myogenic mechanism of each smooth muscle cell we adopt a highly detailed model of subcellular dynamics that has been previously applied on the renal afferent arteriole myocytes [29]. Additionally, to represent accurately conducted responses that have been shown, experimentally and theoretically, to have a significant effect on overall autoregulation [17,24,27,53,54], we incorporate gap junction coupling directly among neighboring smooth muscle cells as well as among smooth muscle cells and a layer of endothelial cells. Lastly, in order to maintain a low computational load, we adopt only a simplified representation of the endothelium similar to existing approaches [7,27,55].

The combined submodels are described in detail below.

2.1.1. Vascular Blood Flow

To obtain an accurate representation of blood flow, i.e., volume flow and pressure profile, we model a segment of the renal vasculature consisting of an afferent and the associated efferent arterioles as shown in Figure 1. Consistent with the renal microanatomy, we represent the two vessels as consecutive straight tubes that extend from the cortical radial artery ($x=0$) to the glomerulus ($x=L_{\mathrm{AA}}$), and from the glomerulus to the entrance of the peritubular capillary network ($x=L_{\mathrm{AA}}+L_{\mathrm{EA}}$), respectively.

Let $Q_{\mathrm{AA}}$ and $Q_{\mathrm{EA}}$ denote the blood flow in the model afferent and efferent arterioles, respectively. In the glomerulus, which is the initial part of the nephron and where filtration takes place [1], a portion of the supplied blood is removed from the bloodstream. To model glomerular filtration, we assume a constant filtration fraction $f_{\mathrm{g}}$. Hence, the single nephron glomerular filtration rate equals $f_{\mathrm{g}}{Q}_{\mathrm{AA}}$ and conservation of mass reads

$$Q_{\mathrm{EA}}=(1-f_{\mathrm{g}})Q_{\mathrm{AA}}.$$

(1)

Blood flow in the model vasculature is assumed to obey Poiseuille’s law [56], which, applied to the two vessels, reads

$$\begin{array}{cc}\hfill \frac{\partial P}{\partial x}& =-\frac{8\mu }{\pi R^4}{Q}_{\mathrm{AA}},0<x<L_{\mathrm{AA}},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial P}{\partial x}& =-\frac{8\mu }{\pi R^4}{Q}_{\mathrm{EA}},L_{\mathrm{AA}}<x<L_{\mathrm{AA}}+L_{\mathrm{EA}}.\hfill \end{array}$$

(3)

In the above equations, P and R denote the pressure and radius profiles along the model vasculature, respectively, and $\mu$ denotes the viscosity of blood, which is assumed equal in both vessels despite a minor difference in the hematocrit caused by filtration [57].

Micro-puncture studies indicate minor pressure differences between the renal artery and the entrance of the afferent arteriole, as well as between the peritubular capillaries and the renal vein [58,59]. Therefore, we use the boundary conditions

$$\begin{array}{cc}\hfill P(t,0)& =P_{\mathrm{a}},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill P(t,L_{\mathrm{AA}}+L_{\mathrm{EA}})& =P_{\mathrm{v}},\hfill \end{array}$$

(5)

where $P_{\mathrm{a}}$ and $P_{\mathrm{v}}$ denote the blood pressures in the renal artery and vein, respectively. As in previous modeling studies [7,55,60], we assume $P_{\mathrm{a}}$ is a prescribed function of time serving as the overall input to our model, while $P_{\mathrm{v}}$ is fixed. To facilitate the presentation in Section 3, we refer to the pressures $P_{\mathrm{AA}}(t,0)$ and $P_{\mathrm{AA}}(t,L_{\mathrm{AA}})$ as afferent arteriole inflow and outflow pressures, respectively.

Conservation of mass (1), Poiseuille’s law (2) and (3), and the boundary conditions (4) and (5) can be combined to yield

$$Q_{\mathrm{AA}}=\frac{P_{\mathrm{a}}-P_{\mathrm{v}}}{W_{\mathrm{AA}}+W_{\mathrm{EA}}},$$

(6)

where the vascular resistances $W_{\mathrm{AA}}$ and $W_{\mathrm{EA}}$ are given by

$$\begin{array}{cc}\hfill W_{\mathrm{AA}}& =\frac{8\mu }{\pi }{\int }_0^{L_{\mathrm{AA}}}\frac{1}{R^4}dx,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill W_{\mathrm{EA}}& =(1-f_{\mathrm{g}})\frac{8\mu }{\pi }{\int }_{L_{\mathrm{AA}}}^{L_{\mathrm{AA}}+L_{\mathrm{EA}}}\frac{1}{R^4}dx.\hfill \end{array}$$

(8)

Autoregulatory phenomena in the renal vasculature of mammals typically require 1–60 s to develop [2,3,4]. On this time scale, the radius of the efferent arteriole appears essentially constant [3], therefore for $x>L_{\mathrm{AA}}$ we assume $R(t,x)=R_{\mathrm{EA}}$, where $R_{\mathrm{EA}}$ is a fixed radius. So (8) reduces to

$$W_{\mathrm{EA}}=\frac{8(1-f_{\mathrm{g}})\mu L_{\mathrm{EA}}}{\pi R_{\mathrm{EA}}^4}.$$

(9)

In contrast, over the same time scale, the radius of the afferent arteriole may change considerably due to either spontaneous contractions of the arteriolar smooth muscles or the operation of the autoregulatory mechanisms [6,7]. Therefore, to accurately predict the evolution of $R(t,x)$ for $x<L_{\mathrm{AA}}$, we develop a detailed model of the afferent arteriolar wall, including the constituting myocytes, which is described below.

2.1.2. Vascular Wall

The afferent arteriole wall model consists of a chain of $N_{\mathrm{AA}}$ smooth muscle cells oriented circumferentially along the vascular lumen as shown in Figure 2. In our formulation, the muscle cells are centered at

$$x_i=\left(i-\frac{1}{2}\right)h,i=1,\dots ,N_{\mathrm{AA}},$$

(10)

where $h=L_{\mathrm{AA}}/N_{\mathrm{AA}}$ is the axial muscle width, i.e., individual smooth muscles are located between $x_i-h/2$ and $x_i+h/2$.

For each smooth muscle cell we adopt the model of [29] which represents in detail subcellular processes such as transmembrane ionic transport, Ca${}^{2+}$ dynamics, and kinetics of myosin light chain phosphorylation. For convenience, we include a schematic of the myogenic mechanism in Figure 3, and refer the reader to [29,30] for complete details on model equations. Among other variables, the smooth muscle cell model predicts muscle membrane potential $v_i$, cytosolic concentrations ${\left[k\right]}_i$ for $k=$ K${}^{+}$, Na${}^{+}$, Cl${}^{-}$, Ca${}^{2+}$, and the fraction ${\psi }_i$ of phosphorylated myosin light chains.

Given our focus on the representation of vascular responses, which are greatly influenced by intercellular communication [17,27,53,55,61,62], we modify the model of [29] which studies cells in isolation and so it ignores signals shared between different smooth muscles. In the present study we consider the addition of gap junctional coupling motivated by the models of [27,30,55].

To represent gap junctions, we model two pathways: a direct one where ions pass between the cytosol of adjacent smooth muscles, and an indirect one where electrical currents pass between the smooth muscles and the endothelium [53,61,63,64]. For the former, we compute the ionic fluxes according to the well established Goldman-Hodgkin-Katz model similar to [63]. For the latter, we assume each smooth muscle is associated with a phenomenological endothelial compartment with local membrane potential $u_i$, and compute the electrical fluxes according to Ohm’s law similar to [27]. We note that on the timescales of the myogenic response that we are interested in, i.e., 1–10 s [3,4,7], the endothelium merely acts as a low-resistance pathway between the membranes of distant cells [3,53,62], and thus we formulate only this aspect. Accordingly, in the present study we ignore vasomodulatory effects mediated by the endothelium such as production and secretion of various vasomodulators such as NO [4,5], which in the mammalian vasculature develop over considerably slower time scales [3,4]. We note that including a full model of the endothelium, for example similar to [63,64,65,66], would considerably add to the computational cost, without providing further autoregulatory insight. Additionally, the experimental observations that exist currently do not suffice for a more elaborate formulation that is specific to the renal microcirculation.

In summary, the modified muscle dynamics are given by

$$\begin{array}{cc}\hfill C_{\mathrm{m}}\frac{d{v}_i}{dt}& =-I_{\mathrm{net}}^i\left(P_i\right)-I_{\mathrm{m}\leftrightarrow \mathrm{m}}^i-I_{\mathrm{m}\leftrightarrow \mathrm{e}}^i\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\mathrm{V}}_{\mathrm{cyt}}\frac{d{\left[k\right]}_i}{dt}& =J_{\mathrm{net},k}^i\left(P_i\right)+J_{\mathrm{m}\leftrightarrow \mathrm{m},k}^i\hfill \end{array}$$

(12)

for $k=$ K${}^{+}$, Na${}^{+}$, Cl${}^{-}$, Ca${}^{2+}$ and the potential dynamics of the endothelial nodes $u_i$ are given by

$$C_{\mathrm{e}}\frac{d{u}_i}{dt}=I_{\mathrm{m}\leftrightarrow \mathrm{e}}^i-I_{\mathrm{e}\leftrightarrow \mathrm{e}}^i$$

(13)

Above, $C_{\mathrm{m}}$ and $C_{\mathrm{e}}$ denote the capacitances of the smooth muscle membrane and the endothelial compartments, respectively; ${\mathrm{V}}_{\mathrm{cyt}}$ denotes the volume of the muscle cytosol; $I_{\mathrm{net}}^i\left(P_i\right)$ denotes the net sum of the currents passing through membrane channels [29]; $I_{\mathrm{m}\leftrightarrow \mathrm{m}}^i$, $I_{\mathrm{e}\leftrightarrow \mathrm{e}}^i$, and $I_{\mathrm{m}\leftrightarrow \mathrm{e}}^i$ denote the currents passing through gap junctions developed between muscles and the endothelium [27,63]; and $J_{\mathrm{net},k}^i$ and $J_{\mathrm{m}\leftrightarrow \mathrm{m},k}^i$ denote the corresponding ion fluxes [29,63]. For the above equations, we assume no-flux boundary conditions at both end points of the afferent arteriole.

The model muscle cells represent the myogenic response under the assumption that the afferent arteriole responds directly to changes in blood pressure [6]. For this, (11) and (12) incorporate in $I_{\mathrm{net}}^i\left(P_i\right)$ and $J_{\mathrm{net},\mathrm{k}}^i\left(P_i\right)$, respectively, contributions from pressure sensitive membrane channels whose activation depends on local blood pressure

$$P_i\left(t\right)=P(t,x_i),i=1,\dots ,N_{\mathrm{AA}}$$

(14)

where $P(t,x)$ is the pressure profile of (2).

We denote the fraction of myosin light chains that are phosphorylated by ${\psi }_i$ and refer to [29,30] for its dependence on the model variables. Given ${\psi }_i$ and pressure $P_i$, the model of [29] predicts the local vascular radius $R_i$ in the vicinity of the ith smooth muscle. In particular, $R_i$ is determined by the balance of apparent hoop stresses $T_P^i$ and $T_{\mathrm{wall}}^i$ (actual hoop stresses are proportional to $T_P^i$ and $T_{\mathrm{wall}}^i$) that are exerted across the muscle,

$$\eta \frac{d{R}_i}{dt}=T_P^i(P_i,R_i)-{\xi }_i{T}_{\mathrm{wall}}^i({\psi }_i,R_i)$$

(15)

The transmural pressure $T_P^i$ depends on $P_i$ and tends to stretch the smooth muscle, causing passive vasodilation. The wall stress $T_{\mathrm{wall}}^i$ depends on ${\psi }_i$ and tends to compress the muscle, causing active vasoconstriction [29]. As the pressure $P_i$ decreases along the vessel [55], the stress $T_P^i$ is lower near the glomerulus than near the cortical radial artery. To achieve a baseline radius profile that is approximately flat, similar to [7,27,60], we introduce a parameter ${\xi }_i$ in (15) which downscales $T_{\mathrm{wall}}^i$ similar to a baseline $T_P^i$. That is, we set

$$\begin{array}{c}\hfill {\xi }_i=1+\left(\frac{P_{\mathrm{g}}^{\mathrm{ref}}}{P_{\mathrm{a}}^{\mathrm{ref}}}-1\right)\frac{x_i}{L_{\mathrm{AA}}}\end{array}$$

(16)

where $P_{\mathrm{a}}^{\mathrm{ref}}$ and $P_{\mathrm{g}}^{\mathrm{ref}}$ denote reference values for the afferent arteriole inflow and outflow pressures.

2.2. Numerical Methods

To obtain solutions of the model equations presented in the previous section we need to solve a large system of (stiff) ordinary differential equations that describes the combined dynamics of the individual smooth muscle cells. This system can be put in the form

$$\frac{d{X}_i}{dt}=F_i(X_{i-1},X_i,X_{i+1},P_i),i=1,\dots ,N_{\mathrm{AA}}$$

(17)

where $X_i$ combines the state variables of the $i^{\mathrm{th}}$ smooth muscle cell (for example $X_i$ combines $v_i$, $u_i$, ${\left[k\right]}_i$, ${\psi }_i$, $R_i$, etc.). We implement no-flux boundary conditions assuming ghost cells $X_0=X_1$ and $X_{N_{\mathrm{AA}}+1}=X_{N_{\mathrm{AA}}}$.

The system in (17) is coupled to the blood flow representation of Section 2.1.1, which is discretized spatially. For this we assume that the radius profile is locally approximated by the smooth muscle predictions

$$R(t,x_i)=R_i\left(t\right),i=1,\dots ,N_{\mathrm{AA}}.$$

(18)

Given (6)–(8), blood flow in the afferent arteriole is computed by

$$Q_{\mathrm{AA}}=\frac{\pi }{8\mu }\frac{P_{\mathrm{a}}-P_{\mathrm{v}}}{{\sum }_i\frac{h}{R_i^4}+(1-f_{\mathrm{g}})\frac{L_{\mathrm{EA}}}{R_{\mathrm{EA}}^4}}.$$

(19)

Pressures at the discrete locations $P_i$ are obtained via (2) using an upwind approximation of the gradient $\partial P/\partial x$. Namely, the pressures are given by

$$\begin{array}{cc}\hfill P_1& =P_{\mathrm{a}}-\frac{4\mu h}{\pi }\frac{Q_{\mathrm{AA}}}{R_1^4},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill P_{i+1}& =P_i-\frac{4\mu h}{\pi }\left(\frac{1}{R_{i-1}^4}+\frac{1}{R_i^4}\right)Q_{\mathrm{AA}},i=2,\dots ,N_{\mathrm{AA}}.\hfill \end{array}$$

(21)

To compute numerical solutions we cast the resulting system of semi-discrete Equations (17) and (19)–(21) in the differential-algebraic form

$$\begin{array}{cc}\hfill \frac{dY}{dt}& =G_1(t,Y,Z),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill 0& =G_2(t,Y,Z),\hfill \end{array}$$

(23)

where $Y=(X_1,\dots ,X_{N_{\mathrm{AA}}})$ and $Z=(Q_{\mathrm{AA}},P_1,\dots ,P_{N_{\mathrm{AA}}})$. For the time evolution of (22) and (23), we apply standard numerical methods for initial value problems in differential-algebraic form [67,68] implemented in MATLAB (The MathWorks, Natick, MA) using implicit solvers that can efficiently treat the stiffness inherent in (17). Specifically, we opt for the solver ode15s, and supply the Jacobian matrix to improve the solver’s efficiency. We note that this built-in MATLAB method is adaptive, so that we expect that the numerical error arising from time discretization is controlled. The absolute and relative tolerances we use for ode15s are both ${10}^{-7}$.

2.3. Parameter Values

2.3.1. Vascular Geometry and Hemodynamics

Throughout this study, we consider an afferent arteriole of total length $L_{\mathrm{AA}}$ = 60 $\mathsf{\mu }$m that consists of cells with axial length h = 3 $\mathsf{\mu }$m approximately of the same dimensions as reported in [69,70]. Based on the estimated $L_{\mathrm{AA}}$ and h, the model afferent arteriole is discretized into $N_{\mathrm{AA}}=20$ cells. We note that while the axial length h used is realistic, in order to keep the overall computational costs low, the vessel length $L_{\mathrm{AA}}$ used is near the smallest measured experimentally in the rat kidney.

Values for the model parameters are listed in

Table 1. Model parameter values.

Description	Parameter	Value	Units	Reference
Reference afferent arteriole inflow pressure	$P_{\mathrm{a}}^{\mathrm{ref}}$	100	mmHg	[55,60]
Reference afferent arteriole outflow pressure	$P_{\mathrm{g}}^{\mathrm{ref}}$	50	mmHg	[55,60]
Reference afferent arteriole radius	$R_{\mathrm{AA}}^{\mathrm{ref}}$	10	$\mathsf{\mu }$m	[29,60]
Reference afferent arteriole flow	$Q_{\mathrm{AA}}^{\mathrm{ref}}$	300	nL/min	[55,60]
Reference glomerular filtration rate	$Q_{\mathrm{g}}^{\mathrm{ref}}$	30	nL/min	[7,60]
Pressure in renal vein	$P_{\mathrm{v}}$	4	mmHg	[7]
Axial length of smooth muscle cell	h	3	$\mathsf{\mu }$m	[69]
Number of smooth muscle cells	$N_{\mathrm{AA}}$	20	-	present study
Afferent arteriole length	$L_{\mathrm{AA}}$	60	$\mathsf{\mu }$m	[69]
Efferent arteriole length	$L_{\mathrm{EA}}$	60	$\mathsf{\mu }$m	present study
Efferent arteriole radius	$R_{\mathrm{EA}}$	11	$\mathsf{\mu }$m	[71]
Filtration fraction	$f_{\mathrm{g}}$	0.1	-	[72]
Blood viscosity	$\mu$	6.68	mmHg ·$\mathsf{\mu }$s	present study
Muscle contraction time constant	$\eta$	1.71	s	[29,30]

. These are chosen based on previous modeling studies or have been computed assuming reference physiological values for the afferent arteriole inflow and outflow pressures $P_{\mathrm{a}}^{\mathrm{ref}}$ and $P_{\mathrm{g}}^{\mathrm{ref}}$, afferent arteriole radius $R_{\mathrm{AA}}^{\mathrm{ref}}$, afferent arteriole blood flow $Q_{\mathrm{AA}}^{\mathrm{ref}}$, and single nephron glomerular filtration rate $Q_{\mathrm{g}}^{\mathrm{ref}}$. We set the radius of the efferent arteriole $R_{\mathrm{EA}}$ to 10% larger than $R_{\mathrm{AA}}^{\mathrm{ref}}$ as reported in [71], and we compute the length of the efferent arteriole $L_{\mathrm{EA}}$ such that the resulting pressure drops along the afferent and efferent arterioles at reference are $P_{\mathrm{a}}^{\mathrm{ref}}-P_{\mathrm{g}}^{\mathrm{ref}}$ and $P_{\mathrm{g}}^{\mathrm{ref}}-P_{\mathrm{v}}$, respectively. The value $\mu$ is adjusted to ensure that, for the dimensions of our system and physiological values of the parameters as listed in Table 1, we reproduce a reference pressure drop of $50$mmHg in Equation (19) in agreement with experimental observations [56].

2.3.2. Electrophysiology

Values for the parameters modeling the currents $I_{\mathrm{net}}^i$ and ion fluxes $J_{\mathrm{net},\mathrm{k}}^i$ of (11) and (12) are adopted from [29,30], while values for the parameters modeling $I_{\mathrm{m}\leftrightarrow \mathrm{m}}^i$ and $J_{\mathrm{m}\leftrightarrow \mathrm{m},\mathrm{k}}^i$ are adopted from [63]. For the currents $I_{\mathrm{m}\leftrightarrow \mathrm{e}}^i$ and $I_{\mathrm{e}\leftrightarrow \mathrm{e}}^i$ modeling endothelial gap junctions, we use Ohmic conductances $g_{\mathrm{me}}$ and $g_{\mathrm{ee}}$, respectively. We compute these values based on the estimated dimensions of the model afferent arteriole, namely

$$\begin{array}{cc}\hfill g_{\mathrm{me}}& =\frac{A_{\mathrm{me}}}{{\rho }_{\mathrm{me}}{\ell }_{\mathrm{me}}},\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill g_{\mathrm{ee}}& =\frac{A_{\mathrm{ee}}}{{\rho }_{\mathrm{ee}}{\ell }_{\mathrm{ee}}},\hfill \end{array}$$

(25)

where $A_{\mathrm{me}}=2\pi h{R}_{\mathrm{AA}}^{\mathrm{ref}}$ is the contact area between the muscles and the endothelium, $A_{\mathrm{me}}=\pi R_{\mathrm{AA}}^{\mathrm{ref}}{\ell }_{me}$ is the endothelial cross-section area, ${\ell }_{\mathrm{me}}$ is half of the afferent arteriole wall thickness, and ${\ell }_{\mathrm{ee}}$ equals h. For the resistivities ${\rho }_{\mathrm{me}}$ and ${\rho }_{\mathrm{ee}}$ we adopt the values from [65].

3. Model Results

Having described the formulation of our model and the necessary numerical methods in the earlier sections, here we show characteristic simulations. Specifically, we examine the model’s behavior when the inflow pressure, denoted $P_{\mathrm{a}}\left(t\right)$, is (i) set to a constant value (time-independent); (ii) perturbed by an instantaneous rise or drop (time-dependent); and (iii) perturbed by sinusoidal oscillations (time-dependent). We focus our presentation only on $P_{\mathrm{a}}\left(t\right)$, which essentially represents the renal perfusion pressure, due to its physiological and experimental importance [2,3,4,5,56].

3.1. Responses to Steady Perturbations

We first compare our model to the previous models of [29] and [30] by studying the base case situation when the inflow pressure $P_{\mathrm{a}}$ is set to a constant value of 100 mmHg (baseline). Figure 4A shows predicted oscillations for the first cell in the afferent arteriole for the membrane potential $v_1$. The mean value is $-35.96$ mV, which is similar to the values found in both [29,30] and which approximates the measured value of $-40$ mV well [73]. Figure 4B shows the cytosolic Ca${}^{2+}$ concentration, which varies between approximately 200 and 360 nM, a somewhat larger range than as predicted in [29,30]. The frequency of the oscillations is 0.15 Hz, the same as in [30] and slightly smaller compared to [29], but well within the range of experimental measurements, for example see [26]. Figure 4C shows the local vascular diameter in the proximity of the first cell of the afferent arteriole, which has an average of 19.7 $\mathsf{\mu }$m with an oscillation amplitude of 0.6 $\mathsf{\mu }$m. Both values are in agreement with previous model predictions [29] and also experimental observations of spontaneous vasomotion, for example [4,71].

To investigate the myogenic response of the model’s afferent arteriole, we computed the outflow pressure and time- and space-averaged diameter of the afferent arteriole given different constant inflow pressures $P_{\mathrm{a}}$. For a comparison, the dotted line in Figure 5 indicates the values expected in the absence of autoregulation (i.e., inflow perturbations are transmitted unattenuated to outflow perturbations) and the dashed line indicates the values expected in the presence of perfect autoregulation (i.e., inflow perturbations are not transmitted to outflow perturbations). We emphasize that, as these lines are shown only for visual purposes, their computation is not based on the present model and instead they can be directly derived from basic principles [56]. As seen in Figure 5, our results indicate that as $P_{\mathrm{a}}$ increases, the predicted outflow pressure of the afferent arteriole increases only when $P_{\mathrm{a}}$ is below ≈80 mmHg or above ≈180 mmHg. The predicted pressure response lies below the line representing perfect autoregulation for low $P_{\mathrm{a}}$, coincides with the line for a wide range of inflow pressures (≈80–180 mmHg), and then lies above the line for large $P_{\mathrm{a}}$. A similar autoregulatory plateau is predicted for $Q_{\mathrm{AA}}$ (not shown). We note that this plateau is in almost complete agreement with Figure 8 in [11], where renal blood flow is measured in vivo over a range of arterial blood pressures in whole kidney preparations. Our model is able to achieve almost perfect autoregulation for a physiologically relevant range of inflow pressures that is larger than the corresponding range predicted by the cellular models of [30,55].

Figure 6A shows the steady-state diameters (time- and space-averaged across afferent arteriole cells) for different constant inflow pressures. As in [30,55], the model predicts vasodilation for small inflow pressures and vasoconstriction for large inflow pressures. Our model predicts a larger mean diameter for small inflow pressures (≤80 mmHg) than [30,55] and approximately the same mean diameter for large inflow pressures (≥160 mmHg).

As can be seen from Figure 6A, vasoresponses to either reduced of elevated inflow pressure are considerably stronger near the distal part of the model afferent arteriole. For example, at 180 mmHg, vasoconstriction of the distal cells is nearly twice that of the proximal cells. For all inflow pressures, according to Figure 6B, attenuation of the perturbations occurs throughout the length of the model afferent arteriole. As a result of the stronger distal responses, the attenuation becomes gradually stronger as perturbations travel downstream.

3.2. Responses to a Step Perturbation

To gain insights into the behavior of the model under time-dependent pressure perturbations we simulate changes in inflow pressure as illustrated in Figure 7, first row. We simulate a pressure pulse by introducing an almost instantaneous rise (respectively drop) in inflow pressure, maintaining the pulse pressure for 20 s, and returning back to the baseline through an almost instantaneous drop (respectively rise). The system is then allowed to stabilize after the return to the baseline inflow pressure over a period of 50 s. This framework allows us to analyze the predicted responses of the model afferent arteriole over a short-term period for both increases and decreases in blood pressure and is intended to simulate step-pressure experimental conditions [74].

Rows 2–5 of Figure 7 show the evolution of flow $Q_{\mathrm{AA}}$, membrane potential, cytosolic Ca${}^{2+}$ concentration, and local vascular diameter for the first cell in the modeled afferent arteriole under a pressure pulse decrease to 80 mmHg (left) and a pressure pulse increase to 120 mmHg (right). The pressure decrease leads to a reduction in membrane potential, which triggers a rapid decrease in the amplitude of the Ca${}^{2+}$ concentration oscillations and subsequently a slower drop in the local vascular diameter [29,30]. This reduction in the local vascular diameter corresponds to passive constriction of the model afferent arteriole due to the sudden pressure drop, which is followed by a slower dilation due to the activation of the myogenic response (Figure 7A5) [29]. The short pressure pulse does not allow the local vascular diameter to fully dilate to a diameter larger than the baseline 20.1 $\mathsf{\mu }$m during the pressure decrease. However, the local vascular diameter returns to a larger mean amplitude of oscillations after the return to baseline pressure at 40 s that stabilizes slowly (≫40 s) to the pre-step value. Qualitatively, our model predicts correctly this active response [3,4]. Furthermore, although somewhat slower than reported in the renal experimental literature [3,4], our predicted time course is also quantitatively similar to the predictions of existing vascular models, for example [23,75,76].

In the pressure pulse increase, pressure-activated Ca${}^{2+}$ and Na${}^{+}$ membrane channels open up, allowing the influx of cations, which in turn lead to an increase in membrane potential as well as cytosolic Ca${}^{2+}$ concentration [29,30]. Following the initial passive dilation, this triggers a rapid active constriction (Figure 7B4, at 20 s). The short pulse period does not allow us to visualize complete constriction in the case of this pressure up-step.

To investigate the propagation of the pressure pulses (Figure 7, first row) throughout the length of the model afferent arteriole, we show the time evolution of membrane potential, cytosolic Ca${}^{2+}$ concentration, and local vascular diameter of 5 different cells located throughout the length of the model vessel during a pressure pulse increase to 140 mmHg in Figure 8. As a result of including the gap junctions (see Section 2.1.2), the membrane potential profiles are similar for all cells. On the other hand, the Ca${}^{2+}$ concentration oscillation amplitude after the pulse shows that the distal cells (cells spanning the terminal ∼20% of the afferent arteriole) differ considerably from their neighbors. The Ca${}^{2+}$ fluctuations further affect the local muscle dynamics such that the diameters of the last model afferent arteriole cells also illustrate larger differences from the proximal and middle cells.

Figure 9 further illustrates the effect of the myogenic response on the outflow diameter and pressure for the short pulse simulations. Figure 9A corresponds to the steady-state diameter of the afferent arteriole (calculated by time- and space- averaging across all cells). These predictions show that some vasodilation at small inflow pressure pulses and some vasoconstriction at large inflow pressure pulses are observed even when investigating short-term responses to pressure changes. Figure 9B shows that a similar behavior is observed for the time-averaged diameter of the last cell in the model afferent arteriole, as predicted in [55] as well. These predictions indicate that the autoregulatory plateau is conserved, though to a lesser degree, even when considering short-term responses.

3.3. Responses to Sinusoidal Perturbation

Motivated by experimental studies reporting influences of oscillatory pressure on blood flow [20], we next examined the model afferent arteriole’s behavior under sinusoidal perturbations in the inflow pressure $P_{\mathrm{a}}\left(t\right)$. After setting the inflow pressure to a constant value of $P_{\mathrm{a}}^{\mathrm{ref}}$ for $t<0$, we perturb the inflow pressure in the form

$$P_{\mathrm{a}}\left(t\right)=P_{\mathrm{a}}^{\mathrm{ref}}+A\sin \left(2\pi ft\right),t\ge 0,$$

(26)

where A is the amplitude of the sinusoidal inflow pressure and f is the frequency. We studied the model’s response to the perturbations for amplitudes $A=10$ mmHg and $A=20$ mmHg and for a range of frequencies between 0.01 Hz (slower) and 1 Hz (faster).

For each inflow amplitude A and frequency f, we observe that the oscillations in the outflow pressure contain two or more different frequencies. Three such examples are given in Figure 10, where the first row shows the perturbed inflow pressure, the second row shows oscillations in the local vascular diameters for the first and last cells in the afferent arteriole, and the third row shows the resulting outflow pressure. To facilitate the comparison between the induced $P_{\mathrm{a}}\left(t\right)$ and predicted $Q_{\mathrm{AA}}\left(t\right)$ perturbations we normalize the corresponding time courses by their respective baseline values and combine them in Figure 11.

The larger the frequency f, the more distinct that a “faster oscillation” with a larger frequency as well as a “slower oscillation” with a smaller frequency appears in the outflow pressure curves (Figure 10, third row). The faster oscillation has approximately the same frequency as the inflow frequency f. As f increases the frequency of the slower oscillation increases and then plateaus around the value of approximately 0.15 Hz for $f\ge 0.5$ Hz (for both amplitudes $A=10$ mmHg and $A=20$ mmHg), which is the natural frequency of the afferent arteriole spontaneous vasomotion reported in Section 3.1.

The amplitude of the outflow pressure is larger or equal to the inflow amplitude A for small frequencies and decreases as the inflow frequency f increases. When the inflow frequency $f\ge 0.25$ Hz, for $A=10$ mmHg, the outflow amplitude converges to approximately $7.8$ mmHg, and for $A=20$ mmHg, the outflow amplitude converges to approximately $13.5$ mmHg. Thus, for smaller inflow frequencies f, the sinusoidal perturbation leads to irregular oscillations in the local vascular diameters and outflow pressure, and for large inflow frequencies f, the perturbation leads to sustained vasoconstriction [6,20,21,55]. According to the analyses in [26], this behavior is expected since, due to absence of tubuloglomerular feedback (TGF) in the present formulation, dynamic autoregulation at low frequencies f is not incorporated.

Assuming that the afferent arteriole behaved like a rigid tube, a linear relationship between $P_{\mathrm{a}}$ and $Q_{\mathrm{AA}}$ is expected. To see this, we note that (7) reduces to a constant resistance $W_{\mathrm{AA}}$, and therefore (6) implies the linear dependence. For example, under the assumption of a rigid afferent arteriole, a 10% increase of inflow pressure $P_{\mathrm{a}}$, as in Figure 10 and Figure 11, would result in a 10% increase of $Q_{\mathrm{AA}}$ and the sinusoidal oscillations in $P_{\mathrm{a}}$ and $Q_{\mathrm{AA}}$ would be in phase. However, the afferent arteriole is not a rigid tube, so our simulations indicate a non-linear relationship between these variables. We find that, similar to the outflow pressure amplitude, the amplitude of the normalized blood flow $Q_{\mathrm{AA}}/Q_{\mathrm{AA}}^{\mathrm{ref}}$ decreases as the inflow frequency f increases and plateaus to approximately $0.17$ for $A=10$ mmHg and approximately $0.29$ for $A=20$ mmHg. Thus, a 10% increase in inflow pressure results in a 17% increase in blood flow rate, and a 20% increase in inflow pressure results in a 29% increase in blood flow rate. The blood flow leads the inflow pressure (Figure 11) since our model does not have a venous system and we have a rigid tube at the end of the afferent arteriole, not a compliant tube. For small input frequencies f, $Q_{\mathrm{AA}}$ has a phase shift in front of $P_{\mathrm{a}}$, but as f increases, $Q_{\mathrm{AA}}$ is shifted after $P_{\mathrm{a}}$.

4. Discussion

We have developed a multiscale mathematical model of the myogenic response of the mammalian kidney. The model spans subcellular, cellular, and vascular aspects and represents detailed Ca${}^{2+}$ trafficking in each of the afferent arteriole myocytes, as well as the kinetics of myosin light chain phosphorylation, the mechanics of the vascular wall, and blood flow. The developed model is an extension of our preliminary arteriolar model [30], which ignores conduction via the endothelial layer. As of that, the model in [30] might underestimate the effects of conducted myogenic responses [53] and the impact on the overall autoregulatory behavior of the model vessel [55].

Similar to our preliminary model [30], the model of the present study is constructed by connecting a large number of afferent arteriole smooth muscle cells. However, unlike [30], inclusion of the endothelium and intercellular coupling allows a realistic representation of ion diffusion and electric conduction along the afferent arteriole. Further, compared with [30], a fluid dynamics model is also added to relate fluid pressure, fluid flow, and vascular resistance so as to allow accurate predictions of all variables required for the activation of the myogenic response. It is also worth noting that our present model also overcomes a numerical limitation on cell size that restricted our approach in [30].

Characteristic simulations show that at physiologic blood pressure, the model predicts spontaneous oscillations (Figure 4) in cytosolic [Ca${}^{2+}$] and in vascular diameter consistent with numerous experimental and theoretical findings, for example [2,3,4,26]. Those oscillations arise from the dynamic exchange of Ca${}^{2+}$ between the cytosol and the sarcoplasmic reticulum, coupled to the stimulation of Ca${}^{2+}$-activated potassium and chloride channels, and the modulation of voltage-activated L-type channels [29,30]. Spontaneous oscillations of the afferent arteriole muscle tone result in oscillations in vascular resistance, fluid pressure, and flow [7].

As in predecessor models such as [7,29], the present model formulates the myogenic response under the assumption that the response is initiated by pressure-activated Ca${}^{2+}$ and Na${}^{+}$ channels found in the cell membrane of the afferent arteriole myocytes [62]. However, the precise molecular cascade through which pressure changes activate ion channels in renal vascular smooth muscle cells has not been characterized physically yet [4,5]. Hence, our formulation of this aspect remains only phenomenological at this point, constructed in such a way that predicted integrated vasoresponses agree with experimental findings relating luminal blood pressure, cell membrane potential, and cytosolic Ca${}^{2+}$ and Na${}^{+}$ concentrations [29,55]. Despite this phenomenological assumption, the developed model afferent arteriole is able to stabilize, to a significant degree, outflow pressure for a range of steady-state inflow pressure spanning 80–180 mmHg (Figure 5).

Additionally, the afferent arteriole model described here reproduces the myogenic response even in the dynamic context. Characteristically, the model maintains an autoregulatory plateau even when challenged with short-term pulses in the inflow pressure (Figure 9). Differences in the distal cells in the model vessel are observed, and may be explained by the fluctuations in Ca${}^{2+}$ due to the pressure pulses. Additionally, sinusoidal perturbations in the inflow pressure are also investigated, and the responses are found to depend on the amplitude and frequency of the perturbation. Since the afferent arteriole is an actively compliant tube, the effect of different amplitude oscillations are found to influence both the increase in the blood flow rate and the phase shift between inflow pressure and flow rate.

Among the limitations of the present study is the high computational cost of modeling large vessels consisting of multiple cells. For this reason, we use a simplified endothelium model similar to [7,27,55], although physiologically, a more complex model would be necessary to predict dynamics at slower time scales. The model endothelium compartment does however act as a low-resistance pathway between the membranes of distant cells, as expected, for the time scales of the myogenic response under consideration.

In reality, renal hemodynamics is a significantly more complicated process than modeled here, with dynamics influenced also by inputs from other processes and organs, which we do not currently incorporate into the model. In addition, our model is fully deterministic, and could benefit from a study of more realistic stochastic pressure perturbations. While our results agree with baseline expectations, our work is a modeling study and so it requires experimental validation. Our model allows for carrying out simulations of realistic scenarios and experiments that may not be easily performed otherwise. For example, direct assessment of the myogenic response would require single afferent arteriole or single nephron studies. This requires access to the kidney’s interior, as well as studying afferent arterioles in isolation; these are challenging experiments, partly because preparation procedures may significantly alter the afferent arteriole’s behavior.

Alongside the myogenic response, another major mechanism contributing in renal hemodynamic control is tubuloglomerular feedback (TGF) [1,56]. Independently of the myogenic response which is triggered by local blood pressure perturbations, tubuloglomerular feedback is triggered by salt reabsorption in the downstream nephron that, in turn, depends on filtration rate. As the triggering signals of the two mechanisms are largely independent but interrelated through blood flow which influences both local pressure and filtration rate, highly non-trivial interactions among the two mechanisms are developed [4,7,26,77]. A potential extension of the present afferent arteriole model would be to include a model of nephron transport (e.g., [78,79,80,81,82]) and tubuloglomerular feedback (e.g., [7,83,84,85]). That would result in an integrative model of renal hemodynamics regulation that can be used for studying the interactions between the myogenic and tubuloglomerular feedback mechanisms in the context of renal autoregulation and for investigating changes in renal hemodynamics in pathophysiological conditions, especially under circumstances (e.g., hypertension and diabetes mellitus) involving complex multilevel responses [5].