TY - JOUR
AU - Moulton, D, E
AB - Abstract In this paper, we consider a lumped-parameter model to predict renal pressures and flow rate during a minimally invasive surgery for kidney stone removal, ureterorenoscopy. A ureteroscope is an endoscope designed to work within the ureter and the kidney and consists of a long shaft containing a narrow, cylindrical pipe, called the working channel. Fluid flows through the working channel into the kidney. A second pipe, the ‘access sheath’, surrounds the shaft of the scope, allowing fluid to flow back out of the urinary system. We modify and extend a previously developed model ( Oratis et al., 2018) through the use of an exponential, instead of linear, constitutive law for kidney compliance and by exploring the effects of variable flow resistance, dependent on the presence of auxiliary ‘working tools’ in the working channel and the cross-sectional shapes of the tools, working channel, scope shaft and access sheath. We motivate the chosen function for kidney compliance and validate the model predictions, with ex vivo experimental data. Although the predicted and measured flow rates agree, we find some disagreement between theory and experiment for kidney pressure. We hypothesize that this is caused by spatial pressure variation in the renal pelvis, i.e. unaccounted for in the lumped-parameter model. We support this hypothesis through numerical simulations of the steady Navier–Stokes equations in a simplified geometry. We also determine the optimal cross-sectional shapes for the scope and access sheath (for fixed areas) to minimize kidney pressure and maximize flow rate. 1. Introduction Solid calcifications within the urinary tract have affected mankind for thousands of years, with the oldest known urinary stone discovered in the bladder of an Egyptian mummy (Shattock, 1905). Left untreated, large stones may lead to infection and eventual loss of kidney function and can, on occasion, be life-threatening (Kum et al., 2016). A common method for the destruction and removal of kidney stones is provided by flexible ureterorenoscopy. The procedure is performed by passing auxiliary working tools through the working channel of a ureteroscope, an endoscope designed to visualize and work within the ureter and the kidney (Fig. 1). Via a miniscule light and a camera on the scope tip, the patient’s ureter and kidney are viewed by a urologist. Fig. 1. Open in new tabDownload slide A photograph of a Boston Scientific ureteroscope (left) with the tip of the scope circled and a zoomed-in schematic provided of the scope tip (right). The scope lies within an access sheath, and a working tool sits inside the working channel. A camera and a light are embedded in the scope wall. Dimensions of the scope shaft and working channel are labelled. The blue regions denote areas through which saline passes; |$Q_{\textrm{in}}$| flows into the kidney and |$Q_{\textrm{out}}$| back out through the access sheath. Fig. 1. Open in new tabDownload slide A photograph of a Boston Scientific ureteroscope (left) with the tip of the scope circled and a zoomed-in schematic provided of the scope tip (right). The scope lies within an access sheath, and a working tool sits inside the working channel. A camera and a light are embedded in the scope wall. Dimensions of the scope shaft and working channel are labelled. The blue regions denote areas through which saline passes; |$Q_{\textrm{in}}$| flows into the kidney and |$Q_{\textrm{out}}$| back out through the access sheath. A method for kidney stone destruction is laser lithotripsy, where a long, thin laser fibre is passed through the ureteroscope and uses shock waves to fragment stones. The resulting particles and debris can then obstruct the field of view for the operating surgeon. To mitigate this, a saline solution flows from a bag above the patient, through a long, thin channel within the scope (the working channel) and out of the scope tip. To exit the body, the fluid flows back in the opposite direction along the outside of the scope shaft. There is only a narrow gap between the scope and the ureteral wall, and if too much saline is delivered, overfilling of the kidney (and hence elevated renal pressures) may occur. High kidney pressures during ureteroscopy have been linked to post-operative complications, such as sepsis (Wilson & Preminger, 1990). Thus, to prevent kidney pressures from rising unduly, many urologists insert the scope through a ureteral access sheath—a hollow cylinder that surrounds the scope, temporarily widening the ureter to allow irrigation fluid to flow more freely back out of the kidney alongside the shaft of the scope (Fig. 1). If more flow is required in a procedural setting to clear the area in front of the camera, the typical practice is to increase the inlet pressure by raising or pressurizing the bag of saline solution in an ad hoc manner. It is clear that flow resistance, both through the working channel and back through the sheath, will play a crucial role in determining whether kidney pressures will remain below a physically safe threshold. Resistance to flow can be altered through modifications to the geometry of the surgical tools; in particular, the cross-sectional shapes of the working tools, channel, scope shaft and ureteral access sheaths (Williams et al., 2020). Hence, understanding the relationships between the size and shape of ureteroscopic instruments, inlet pressure, flow rate and pressure within the kidney, has enormous potential to enhance surgical efficiency and safety as well as to revolutionize ureteroscope design. Recent efforts have examined the effects of geometric parameters of the ureteroscope and working tools on components of irrigation flow with a mathematical modelling approach (Williams et al., 2019b). Flow through an isolated ureteroscope has been shown, through asymptotic analysis on the small aspect ratio of the working channel, and comparison with bench-top experiments, to obey a Poiseuille law for flow through a cylindrical pipe. Similarly, flow through a working channel with a working tool, or back through an access sheath, can be modelled as flow through a pipe containing a solid, coaxial cylindrical rod. Flow for a fixed pressure drop through these configurations is lowest when the solid rod is in the centre and increases monotonically with rod proximity to the edge of the bounding channel. This highlights the importance of cross-sectional geometry on viscous resistance, motivating the discussion of optimally shaped elliptical working tool/channel (or similarly scope shaft/access sheath) configurations, which have been theorized to further decrease viscous resistance when compared with circles of equal cross-sectional area (Williams et al., 20120). In this article, we will extend the previous work on component modelling of ureteroscope irrigation and the effects of viscous resistance, to consider a lumped-parameter system model, connecting flow through the working channel and back through the access sheath. Through this, we can understand how geometric changes to particular components of the irrigation system will affect flow rate and kidney pressure. Fluid is driven through the working channel via a prescribed inlet pressure, and by modelling the kidney as a compliant body, and equating fluxes in through the scope and back out through the sheath, a lumped-parameter model for flow rate and kidney pressure as functions of time can be derived (Oratis et al., 2018). The resulting ordinary differential equation depends on the ratio between the viscous resistances through the working channel and back through the sheath and the relationship between changes in kidney pressure and volume (compliance). In Oratis et al. (2018), flow through the working channel and access sheath were modelled as Poiseuille flow through an unobstructed pipe and a concentric annular pipe of circular cross-sections, respectively. The importance of low outflow resistance to maintain low intrarenal pressures was discussed in Oratis et al. (2018) only within the context of increasing access sheath size or intermittently withdrawing the scope. We will extend this work by varying the inflow and outflow resistances to represent the addition of working tools, a non-concentric scope shaft, and modifying the cross-sectional shapes of the scope shaft and access sheath. Additionally, in Oratis et al. (2018), the kidney is modelled as linearly elastic with constant renal pelvis stiffness, |$K$|⁠. However, biological tissue is typically shown to be poorly compliant, represented by an increase in the slope of the stress-strain curve at higher strains (Chagnon et al., 2018; Holzapfel, 2001). This is often accounted for by constitutive laws that include exponential terms, ensuring that tissue stiffens as it stretches. An experimental study (Cruces et al., 2014) performed on porcine kidneys, measuring pressure within the parenchyma in response to fluid volume injected into the kidney, also suggests an exponential relationship between kidney pressure and volume. We will validate the use of an exponential function for the kidney compliance with the use of ex vivo experimental data and will consider its effect on the predictions of the theoretical model. In Section 2, we describe the ex vivo experimental setup and procedure. In Section 3, we derive a lumped-parameter model for kidney pressure and flow rate, considering an exponential renal stiffness, and present the analytical solution. The results in Section 4 are split into three main components. Firstly, in Section 4.1, we validate the use of an exponential renal compliance by fitting this constitutive law to the experimental data and discuss its implications. In Section 4.2, we test the steady-state model predictions for flow rate and kidney pressure against the experimental data. This motivates a discussion of the spatial variance of kidney pressure that is not accounted for in the lumped-parameter model. Finally, in Section 4.3, we theorize on the addition of working tools in the working channel and the effects of changing the geometries of the access sheath and scope shaft on flow rate and kidney pressure. 2. Ex vivo experiments and data extraction A suite of ex vivo experiments was performed on porcine kidneys, with a laboratory configuration shown in Fig. 2(a). The experiments were conducted at room temperature and within several hours of harvesting the kidneys. Irrigation fluid flowed from a saline bag 60 cm above the height of the kidney, when measured from the base of the bag. A pressure cuff surrounded the bag, allowing the inlet pressure to be increased above hydrostatic, and a pressure probe was placed at the inlet to the scope. A 36-cm long 11/13 F Navigator access sheath was inserted through the ureter, with a pressure probe passed between the ureter and the wall of the sheath. A LithoVue ureteroscope was placed within the sheath, and the tips of the sheath, scope and pressure probe were all within the renal pelvis. The kidney and a dish to collect the outflow from the sheath were on OHAUS Traveler mass balances, and FISO Fiber Optic Pressure Sensors recorded the inlet and intrarenal pressures. The pressure sensors were calibrated to atmospheric pressure. Fig. 2. Open in new tabDownload slide Figure (a) shows a photograph of the experimental setup, with measured variables labelled, along with the lengths of the scope and the access sheath. Figures (b) and (c) show examples of data collected from kidney |$A$| for data sets I and II, respectively. Fig. 2. Open in new tabDownload slide Figure (a) shows a photograph of the experimental setup, with measured variables labelled, along with the lengths of the scope and the access sheath. Figures (b) and (c) show examples of data collected from kidney |$A$| for data sets I and II, respectively. Inflations of the pressure cuff surrounding the irrigation bag enabled control of the inlet pressure, and the resulting inlet pressure, kidney mass, kidney pressure and outflow mass were all measured over time. Pressure data were collected at 8-ms intervals, and mass data collected at approximately 500-ms intervals. Experiments were performed on three different kidneys, which we refer to as kidneys A, B and C. Multiple experiments were performed in succession—without adjusting the positions of the experiment equipment—on each of the three kidneys. The kidneys were drained prior to the experiments, so the initial renal pressure before the insertion of irrigation fluid was approximately atmospheric. The inlet pressure was increased over time in two ways: (I) continuously and (II) step-wise, maintaining the pressure of the cuff for approximately 30 seconds before further inflation. Examples of the resulting inlet pressure over time for data sets I and II are shown in Fig. 2(b) and 2(c), respectively. A number of data sets I and II collected from each kidney are summarized in Table 1, and more details on the data analysis are provided in Appendix A. Table 1 Number of experiment runs on each kidney, for each type of data set . A . B . C . Data set I 1 0 3 Data set II 4 2 1 . A . B . C . Data set I 1 0 3 Data set II 4 2 1 Open in new tab Table 1 Number of experiment runs on each kidney, for each type of data set . A . B . C . Data set I 1 0 3 Data set II 4 2 1 . A . B . C . Data set I 1 0 3 Data set II 4 2 1 Open in new tab The two data sets enable us to analyse two distinctly different features: (i) the rise of pressure to steady state and (ii) flow and pressure conditions at steady state for varying inlet pressure. Since data set I consisted of a continuous rise in pressure, this data set was well suited for evaluating the pressure rise and estimating kidney compliance, by extracting a continuous curve (parametrized by time) relating fluid mass and pressure within the kidney. In data set II, since the inlet pressure was held constant for several seconds at each incremental value, this enabled the system to reach steady-state flow (flow into the kidney equalling flow out of the kidney). Thus, this data set was used to study the steady-state conditions, by extracting 5 seconds of data from each interval of constant inlet pressure and using this to determine the steady-state kidney pressure and flow rate as functions of inlet pressure.1 It is also important to note that the two regimes, (i) and (ii), are fundamentally separated in terms of the modelling. In particular, the steady-state conditions are independent of the compliance, as we shall show below. Therefore, the compliance parameter fitting related to (i) is irrelevant to the analysis of (ii). Details on calculating error bars for all data figures are provided in Section A.4. 3. Theoretical formulation We consider irrigation flow in an idealized setting, replicating the laboratory setup in Fig. 2(a). A scope of length |$L_a$| lies horizontally, with fluid of density |$\rho $| and viscosity |$\mu $| entering the working channel of radius |$a$| at fixed pressure |$P_{\textrm{in}}$|⁠. The tip of the scope is inside the kidney which contains irrigation fluid of mass |$M_{\textrm{k}}$|⁠, volume |$V_{\textrm{k}}$| and pressure |$P_{\textrm{k}}$|—all functions of time. An access sheath of radius |$b$| and length |$L_b$| surrounds the scope, with one end within the kidney and the other exposed to atmospheric pressure, which we take to be zero. We term the fluid volumetric flux into the kidney as |$Q_{\textrm{in}}$| and back out through the sheath as |$Q_{\textrm{out}}$|⁠. 3.1 Governing equations Conservation of mass provides $$\begin{equation} \frac{\textrm{d}M_k}{\textrm{d}t} = \rho(Q_{\textrm{in}} - Q_{\textrm{out}}). \end{equation}$$(3.1) The irrigation fluid is taken to be incompressible, and hence increases in fluid mass within the kidney and the volume of fluid within the kidney are correlated directly $$\begin{equation} M_{\textrm{k}} = \rho V_{\textrm{k}}. \end{equation}$$(3.2) Thus, to model kidney pressure as a function of time with (3.1), it remains to determine the kidney compliance: a constitutive relationship between changes in pressure and changes in volume. Defining |$P_{\textrm{k}}$| as the change in kidney pressure from atmospheric, we take $$\begin{equation} P_{\textrm{k}} = P^{\star}(\exp(C V_{\textrm{k}})-1), \end{equation}$$(3.3) where |$C$| is a constant with units of inverse volume and |$P^{\star }$| is a characteristic pressure scaling. Both |$C$| and |$P^{\star }$| will be determined through a fit to experimental data in Section 4.1. The use of (3.3) is motivated by Cruces et al. (2014), Holzapfel (2001) and Chagnon et al. (2018) and has also been seen to provide a good approximation for bladder filling (Griffiths, 1980). Differentiating (3.3) with respect to time, we obtain $$\begin{equation} \frac{\textrm{d}P_{\textrm{k}}}{\textrm{d}t} = P^{\star}C\exp(CV_{\textrm{k}})\frac{\textrm{d}V_{\textrm{k}}}{\textrm{d}t}, \end{equation}$$(3.4) and using (3.2) and (3.1), we obtain the following non-linear differential equation: $$\begin{equation} \frac{\textrm{d}P_{\textrm{k}}}{\textrm{d}t} = C(P^{\star}+P_{\textrm{k}})(Q_{\textrm{in}}-Q_{\textrm{out}}). \end{equation}$$(3.5) As an initial condition for (3.5), we take $$\begin{equation} P_{\textrm{k}}(0) = 0. \end{equation}$$(3.6) 3.2 Viscous resistance We consider steady, unidirectional flows through the working channel and access sheath. The assumption of unidirectional flow through the working channel has been previously validated through comparison with bench-top experiments (Williams et al., 2019a). Steady flow can be justified from the full Navier–Stokes equations and incompressibility condition by considering flow through a working channel (with a similar discussion for flow through the access sheath). Defining the radius-to-length ratio of the channel as |$\epsilon = a/L_a$| and the Reynolds number of the flow as |$\textrm{Re} = QL_a\rho /\mu \pi a^2$| for a typical flow rate |$Q= \mathcal{O}(10^{0})\textrm{ cm}^3\text{/s}$|⁠, the key parameters governing the behaviour of the fluid are |$\epsilon ^2\textrm{Re}$| and |$a^2\rho /\mu T$|⁠, where |$T$| is a chosen time scale. Here, the relevant time scale for (3.5) is the rise time for kidney pressure (or flow rate) to adjust to steady state after changing the inlet pressure. This can be estimated from the ex vivo data to be |$T = \mathcal{O}(10^1\textrm{ s})$|⁠. With values for |$a$|⁠, |$L_a$|⁠, |$\rho $| and |$\mu $| taken from Table 2, |$\epsilon ^2 = \mathcal{O}(10^{-6})$|⁠, |$\epsilon ^2\textrm{Re} = \mathcal{O}(10^{-1})$| and |$a^2\rho /\mu T = \mathcal{O}(10^{-2})$|⁠. Thus, inertial and time-dependent effects within the fluid can be neglected, and we model flow through the working channel and back through the access sheath as steady, unidirectional flow through pipes of arbitrary cross-sections. Table 2 Experimental parameters . Symbol . Value . Unit . Working channel radius |$a$| |$6.00\times 10^{-2}$| cm Scope shaft radius — |$1.58\times 10^{-1}$| cm Access sheath radius |$b$| |$1.93\times 10^{-1}$| cm Working channel length |$L_a$| |$7.90 \times 10^{1}$| cm Access sheath length |$L_b$| |$3.60 \times 10^{1}$| cm Fluid density |$\rho $| |$1.00 \times 10^{0}$| g/cm|$^3$| Fluid viscosity |$\mu $| |$1.00 \times 10^{-2}$| g/cm s . Symbol . Value . Unit . Working channel radius |$a$| |$6.00\times 10^{-2}$| cm Scope shaft radius — |$1.58\times 10^{-1}$| cm Access sheath radius |$b$| |$1.93\times 10^{-1}$| cm Working channel length |$L_a$| |$7.90 \times 10^{1}$| cm Access sheath length |$L_b$| |$3.60 \times 10^{1}$| cm Fluid density |$\rho $| |$1.00 \times 10^{0}$| g/cm|$^3$| Fluid viscosity |$\mu $| |$1.00 \times 10^{-2}$| g/cm s Open in new tab Table 2 Experimental parameters . Symbol . Value . Unit . Working channel radius |$a$| |$6.00\times 10^{-2}$| cm Scope shaft radius — |$1.58\times 10^{-1}$| cm Access sheath radius |$b$| |$1.93\times 10^{-1}$| cm Working channel length |$L_a$| |$7.90 \times 10^{1}$| cm Access sheath length |$L_b$| |$3.60 \times 10^{1}$| cm Fluid density |$\rho $| |$1.00 \times 10^{0}$| g/cm|$^3$| Fluid viscosity |$\mu $| |$1.00 \times 10^{-2}$| g/cm s . Symbol . Value . Unit . Working channel radius |$a$| |$6.00\times 10^{-2}$| cm Scope shaft radius — |$1.58\times 10^{-1}$| cm Access sheath radius |$b$| |$1.93\times 10^{-1}$| cm Working channel length |$L_a$| |$7.90 \times 10^{1}$| cm Access sheath length |$L_b$| |$3.60 \times 10^{1}$| cm Fluid density |$\rho $| |$1.00 \times 10^{0}$| g/cm|$^3$| Fluid viscosity |$\mu $| |$1.00 \times 10^{-2}$| g/cm s Open in new tab These are determined as solutions to $$\begin{equation} \nabla^2w = \frac{1}{\mu}\frac{\textrm{d}p}{\textrm{d}z}, \end{equation}$$(3.7) where |$w$| is the axial velocity as a function of cross-sectional coordinates, and |$\textrm{d}p/\textrm{d}z$| is a constant pressure gradient along the length of the pipe (Lamb, 1895). Flow rate is obtained by integrating the solution to (3.7) for |$w$| over an arbitrary cross-section of the pipe to obtain $$\begin{equation} Q = \frac{1}{\mu}\frac{\textrm{d}p}{\textrm{d}z}\frac{1}{R}, \end{equation}$$(3.8) where |$R$| characterizes the viscous resistance to the geometry of the pipe cross-section. Modelling the working channel as a hollow cylinder of radius |$a$|⁠, $$\begin{equation} Q_{\textrm{in}} = \frac{1}{\mu}\left(\frac{P_{\textrm{in}}-P_{\textrm{k}}}{L_a}\right)\frac{1}{R_{\textrm{in}}}, \qquad R_{\textrm{in}} = 8/\pi a^4. \end{equation}$$(3.9) Assuming that the scope and the sheath are coaxial (although not necessarily concentric), flow back through the access sheath is $$\begin{equation} Q_{\textrm{out}} = \frac{1}{\mu}\frac{P_{\textrm{k}}}{L_b}\frac{1}{R_{\textrm{out}}}, \end{equation}$$(3.10) where |$R_{\textrm{out}}$| can be obtained for circular cylinders as a function of scope position from the analytical solution for flow through an offset annular domain (Heyda, 1959). |$R_{\textrm{out}}$| decreases monotonically as the scope is placed closer to the edge of the sheath. In practice, the position of the scope shaft within the sheath is unknown, so we allow |$R_{\textrm{out}}$| to take any value between its minimum (at the edge of the sheath) and its maximum (in the centre of the sheath). 3.3 Dimensionless model Combining (3.5), (3.9) and (3.10), we obtain $$\begin{equation} \frac{\textrm{d}P_{\textrm{k}}}{\textrm{d}t} = (C/\mu)(P^{\star}+P_{\textrm{k}})\left(\frac{P_{\textrm{in}}-P_{\textrm{k}}}{L_a R_{\textrm{in}}} - \frac{P_{\textrm{k}}}{L_b R_{\textrm{out}}}\right). \end{equation}$$(3.11) We now make the following non-dimensionalizations $$\begin{equation} P_{\textrm{k}} = P^{\star}\hat{P_{\textrm{k}}}, \qquad t = (L_aR_{\textrm{in}}\mu/CP^{\star})\hat{t}, \end{equation}$$(3.12) where hats denote dimensionless variables. The dimensionless form of (3.11) is then $$\begin{equation} \frac{\textrm{d}\hat{P}_{\textrm{k}}}{\textrm{d}\hat{t}} = (1+\hat{P}_{\textrm{k}})(\hat{P}_{\textrm{in}}-(1+\mathcal{R})\hat{P}_{\textrm{k}}), \end{equation}$$(3.13) where we define the dimensionless parameters $$\begin{equation} \hat{P}_{\textrm{in}} = P_{\textrm{in}}/P^{\star}, \qquad \mathcal{R} = (L_aR_{\textrm{in}})/(L_bR_{\textrm{out}}). \end{equation}$$(3.14)|$\hat{P}_{\textrm{in}}$| is the dimensionless inlet pressure and |$\mathcal{R}$| is the ratio of resistance in through the working channel to resistance out through the access sheath. The dimensionless initial condition, (3.6) is $$\begin{equation} \hat{P}_{\textrm{k}}(0) = 0. \end{equation}$$(3.15) 3.4 Analytical solution Equation (3.13) with initial condition (3.15) can be solved by separating variables to obtain the analytical solution $$\begin{equation} \hat{P}_{\textrm{k}} = \frac{(e^{\tau\hat{t}}-1)\hat{P}_{\textrm{in}}}{\hat{P}_{\textrm{in}}+e^{\tau\hat{t}}(1+\mathcal{R})}, \end{equation}$$(3.16) where |$\tau = 1+\hat{P}_{\textrm{in}}+\mathcal{R}$|⁠. The steady-state intrarenal pressure can be easily found by setting the time derivative in (3.13) to zero. We note the presence of |$(1+\hat{P}_{\textrm{k}})$| in (3.13) is due to the chosen exponential function for kidney compliance, |$ \textrm{d}P_{\textrm{k}}/\textrm{d}V_{\textrm{k}}=P_{\textrm{k}}+P^{\star }$|⁠, (3.3). Hence, we can exclude the non-physical steady-state solution |$\hat{P}_{\textrm{k}} = -1$| which corresponds to zero compliance. The relevant steady-state solution is thus $$\begin{equation} \hat{P}_{\textrm{k}}^{\textrm{SS}} = \frac{\hat{P}_{\textrm{in}}}{1 + \mathcal{R}}. \end{equation}$$(3.17) If |$\mathcal{R} \ll 1$|⁠, and hence resistance through the sheath dominates, we see from (3.17) that |$\hat{P}_{\textrm{SS}} \approx \hat{P}_{\textrm{in}}$|⁠. Conversely, if |$\mathcal{R} \gg 1$|⁠, and there is negligible resistance through the sheath, |$\hat{P}_{\textrm{SS}} \approx 0$|⁠. We note that the steady state is not affected by the inclusion of a non-linear kidney compliance, (3.3), and thus the parameters |$C$| and |$P^{\star }$| determine only the time required to reach equilibrium pressure. The steady state will be reached when |$Q_{\textrm{in}} = Q_{\textrm{out}} = Q^{\textrm{SS}}$|⁠, and the dimensional value of the steady-state flow rate can be immediately determined from (3.17) and one of either (3.9) or (3.10) after redimensionalizing and is $$\begin{equation} Q^{\textrm{SS}} = \frac{1}{\mu}\left(\frac{P_{\textrm{in}}}{L_aR_{\textrm{in}}+L_bR_{\textrm{out}}}\right). \end{equation}$$(3.18) 4. Results The parameters used to validate the model are taken from the experimental setup and are summarized in Table 2. 4.1 Estimating compliance (data set I) We validate the use of an exponential compliance (3.3) in Fig. 3(a), by demonstrating a superior fit with experimental data (solid red line) when compared with a linear fit (dashed green line). Here, two data sets were chosen from kidney C, with the rationale behind their use detailed in Appendix A.2; we also note that using kidney A for the data fit results in very similar conclusions (see Fig. A9 in Appendix A.2). The exponential and linear relations are fitted to the mean of the experimental data. From the exponential fit, values are determined as $$\begin{equation} P^{\star} = 5.24 \textrm{ cm H}_2\textrm{O}, \qquad C = 0.43 \textrm{ mL}^{-1}. \end{equation}$$(4.1) Fig. 3. Open in new tabDownload slide (a) Data from kidney C runs (i) and (j) (see Appendix A). Exponential fit (red) with |$P^{\star } = 5.24$| cm H|$_2$|O and |$C = 0.43$| mL|$^{-1}$|⁠. Linear model from Oratis et al. (2018), (dashed blue), with |$K = 0.41$| cm H|$_2$|O/g. Linear fit (dashed green) with |$K=9.59$| cm |$H_20\ \textrm{mL}^{-1}$|⁠. (b) Comparison of theoretical kidney pressure |$P_{\textrm{k}}$| for the three different compliance models from (a). The different time scales required to reach steady state are indicated by the corresponding coloured labels below the |$x$|-axis. The data in Fig. 3(b) is for |$P_{\textrm{in}} = 60$| cm H|$_2$|O. Fig. 3. Open in new tabDownload slide (a) Data from kidney C runs (i) and (j) (see Appendix A). Exponential fit (red) with |$P^{\star } = 5.24$| cm H|$_2$|O and |$C = 0.43$| mL|$^{-1}$|⁠. Linear model from Oratis et al. (2018), (dashed blue), with |$K = 0.41$| cm H|$_2$|O/g. Linear fit (dashed green) with |$K=9.59$| cm |$H_20\ \textrm{mL}^{-1}$|⁠. (b) Comparison of theoretical kidney pressure |$P_{\textrm{k}}$| for the three different compliance models from (a). The different time scales required to reach steady state are indicated by the corresponding coloured labels below the |$x$|-axis. The data in Fig. 3(b) is for |$P_{\textrm{in}} = 60$| cm H|$_2$|O. In Oratis et al. (2018), a linear renal stiffness was assumed with mean value |$K = 0.41$| cm H|$_2$|O mL|$^{-1}$|⁠. This is plotted as the dashed blue line in Fig. 3(a). The shallow slope models the kidney as significantly more compliant than the data suggests; in fact, a linear fit to our data (dashed green line in Fig. 3(a)) produces a renal stiffness values of |$K = 9.59$| cm H|$_2$|O mL|$^{-1}$|⁠. For the red line in Fig. 3(b), we consider a fixed inlet pressure of |$P_{\textrm{in}} = 60$| cm H|$_2$|O and plot (3.16) using the best fit values of |$P^{\star }$| and |$C$|⁠, (4.1). The value for |$R_{\textrm{out}}$| is taken for a scope sitting concentrically in the access sheath. Using the same values for |$P_{\textrm{in}}$|⁠, |$R_{\textrm{in}}$| and |$R_{\textrm{out}}$|⁠, we plot the solution from Oratis et al. (2018) as the dashed blue line. We note that not only do the curves have slightly different shapes but also the time scales required to reach steady state differ significantly. These are indicated by the corresponding coloured labels below the |$x$|-axis. This result is intuitive from the two different compliance relationships, as a stiffer kidney will reach steady state more quickly. Oratis et al. (2018) predicted a characteristic rise time of the order of several minutes, which does not correlate with experimental and clinical observations of |$\mathcal{O}(10^1\textrm{ s})$|⁠. Using the stiffness value from the linear fit (dashed green line in Fig. 3), we get the corresponding dashed green line in Fig. 3(b) for kidney pressure over time. We see that the characteristic rise time is better captured by the fitted linear stiffness value compared with the dashed blue, although the shape of the curve still differs from the red curve corresponding to the exponential fit. In fact, as we are considering kidney pressure |$P_{\textrm{k}} < 30$|⁠, the fitted linear stiffness model estimates the material as stiffer than the more accurate exponential compliance, and thus the rise to steady state is slightly shorter. For higher kidney pressures, the difference in slope will be more marked and the discrepancy between predicted |$P_{\textrm{k}}$| curves even more significant. It is also evident in Fig. 3(a) that the exponential compliance fits the data much more closely, as a linear law cannot capture stiffening effects. 4.2 Testing model predictions (data set II) We now test the steady-state model predictions for kidney pressure and flow rate, (3.17) and (3.18), against the experimental data obtained from data set II. We note that there is no fitting required between the model and the data, as the steady state predictions are independent of kidney compliance. In Fig. 4(a), we plot the steady-state flow prediction, (3.18), along with experimental data, as functions of |$P_{\textrm{in}}$|⁠. We consider both the minimum and maximum values for |$R_{\textrm{out}}$|⁠, when the scope is fully offset within the sheath and centred, respectively (solid black lines), and shade the area between the two curves, corresponding to the predicted |$Q^{\textrm{SS}}$| over the range of feasible |$R_{\textrm{out}}$| values. The experimental flow-rate data for kidneys A, B and C, estimated from the outflow mass values over time, are shown to lie within the predicted range (Fig. 4(a)). Fig. 4. Open in new tabDownload slide (a) Equation (3.18) as the black lines compared with data collected from three different kidneys labelled A, B and C. The dashed red line is the theoretical prediction for the optimal sheath/scope geometry (see Section 4.3); (b) (3.17). The resistance ratio corresponding to a centred scope is |$\mathcal{R} = 1.7$| and to an offset scope is |$\mathcal{R} = 4.1$|⁠. Fig. 4. Open in new tabDownload slide (a) Equation (3.18) as the black lines compared with data collected from three different kidneys labelled A, B and C. The dashed red line is the theoretical prediction for the optimal sheath/scope geometry (see Section 4.3); (b) (3.17). The resistance ratio corresponding to a centred scope is |$\mathcal{R} = 1.7$| and to an offset scope is |$\mathcal{R} = 4.1$|⁠. Accurately predicting the measured steady-state intrarenal pressure with (3.17) is more challenging, and in Fig. 4(b), we plot the wedge of predicted |$P_{\textrm{k}}^{\textrm{SS}}$| and observe that the data from kidney A do not fall within the estimated range. One possible cause of the discrepancy is the position of the pressure probe within the kidney in the experiments; although it was placed proximal to the end of the scope, its position was likely to vary. Within the current modelling framework, we assume pressure within the kidney to be constant, i.e the pressures at the exit of the working channel and entrance to the access sheath are identical and equal to |$P_{\textrm{k}}$|⁠. In actuality, there will be complex flow patterns and pressure variation within the renal pelvis, dependent on flow rate, renal structure and the ureteroscopic instruments in use. To test whether large variations in renal pressure are expected, and whether these can account for the observed discrepancy in Fig. 4(b), we solve the steady Navier–Stokes equations in the simplest representative geometry (see Appendix C for more details). We model the renal pelvis as a 2D rectangular cavity, prescribing a Poiseuille inlet and zero-stress outlet boundary conditions. Defining the Reynolds number as |$\textrm{Re} = \rho Qa/(\mathcal{A}\mu )$|⁠, where |$\mathcal{A} = \pi a^2$| is the cross-sectional area of the working channel, the relevant Reynolds numbers for simulations representing the flow rates in Fig. 4(a), and thus measured renal pressures in Fig. 4(b), range from |$\textrm{Re} = 0$| to |$\textrm{Re} = 1883$|⁠. Results of the simulation for a representative flux |$Q = 1.4$| cm|$^3$|/s (⁠|$\textrm{Re} = 1500$|⁠) are shown in Fig. 5(a) and (b). We consider the rectangular domain as a cross-sectional slice through the renal pelvis; thus, the half-width is taken to be the radius of the access sheath and the radii of the sheath, scope and working channel are taken from the experiments (Table 2) and indicated along the vertical axis of Fig. 5(a) in cm units. The velocity streamlines in Fig. 5(a) demonstrate the presence of multiple vortical structures, with the largest region of closed streamlines near the centre of the cavity. The corresponding pressure contour plot is shown in Fig. 5(b), and the colourbar indicates pressure in units of cm H|$_2$|O. As flow is driven by pressure gradients, the same velocity solution will exist if an arbitrary constant is added to the pressure field; thus, we shift the pressure solution so the minimum value is at zero. Hence, Fig. 5(b) indicates the pressure variation in the cavity, rather than absolute pressure values. We see that the pressure varies significantly on a spatial scale with a maximum variation of |$90$| cm H|$_2$|O. Fig. 5. Open in new tabDownload slide Plot of predicted (a) streamlines and (b) pressure variation, within a rectangular cavity with prescribed parabolic inflow for |$Q = 1.4$| cm|$^3$|/s (⁠|$\textrm{Re} = 1500$|⁠) and zero-stress outflow conditions. Cavity measurements labelled in (a) in cm units, corresponding to working channel, scope and sheath dimensions. Pressure variation in (b) in cm H|$_2$|O, shifted so the minimum pressure in the cavity is at zero. Figure (c) plots the difference between the minimum and maximum pressure values in the cavity, as a function of |$P_{\textrm{in}}$|⁠. The relationship between |$P_{\textrm{in}}$| and |$Q$| prescribed for each simulation is determined from (3.18) with |$R_{\textrm{out}}$| taken as the value for a scope centred in the sheath. Fig. 5. Open in new tabDownload slide Plot of predicted (a) streamlines and (b) pressure variation, within a rectangular cavity with prescribed parabolic inflow for |$Q = 1.4$| cm|$^3$|/s (⁠|$\textrm{Re} = 1500$|⁠) and zero-stress outflow conditions. Cavity measurements labelled in (a) in cm units, corresponding to working channel, scope and sheath dimensions. Pressure variation in (b) in cm H|$_2$|O, shifted so the minimum pressure in the cavity is at zero. Figure (c) plots the difference between the minimum and maximum pressure values in the cavity, as a function of |$P_{\textrm{in}}$|⁠. The relationship between |$P_{\textrm{in}}$| and |$Q$| prescribed for each simulation is determined from (3.18) with |$R_{\textrm{out}}$| taken as the value for a scope centred in the sheath. In Fig. 4(b), the deviation between |$P_{\textrm{k}}$| for kidney A and the predicted values increases with |$P_{\textrm{in}}$|⁠. To estimate |$P_{\textrm{in}}$| for a given simulation with a prescribed flow rate, we use the relationship for the steady-state flow rate, (3.18), assuming |$R_{\textrm{out}}$| for a scope centred within the access sheath. Using this approximation, the simulation in Fig. 4(a) and (b) corresponds to |$P_{\textrm{in}} = 276$| cm H|$_2$|O. We can then plot the maximum pressure variation (difference between maximum and minimum pressures in the cavity) as a function of |$P_{\textrm{in}}$| (Fig. 5(c)), and we see that the variation increases as a function of |$P_{\textrm{in}}$|⁠. In a sense, this predicted variation corresponds to large error bars on the |$P^{\textrm{SS}}_{\textrm{k}}$| data points in Fig. 4(b), explaining our inability to always accurately predict the pressure detected by the probe. Additionally, the pressure variation in the cavity increases with |$P_{\textrm{in}}$|⁠, possibly explaining the increase in deviation between the experiments and prediction for kidney A in Fig. 4(b). 4.3 Effects of working tools and access sheath geometry The presence of working tools, such as laser fibres or stone baskets within the working channel, can be included in (3.17) and (3.18) by increasing |$R_{\textrm{in}}$| according to the size, geometry and placement of the working tool within the working channel. Hence, for fixed |$P_{\textrm{in}}$|⁠, adding a working tool to the channel will decrease |$P^{\textrm{SS}}_{\textrm{k}}$| and hence |$Q^{\textrm{SS}}$|⁠. Kidney pressure as a function of working tool radius (where we recall that the working channel has radius |$0.06$| cm) is plotted in Fig. 6(a) as the dashed black line. As a minimum flow rate is often required to maintain a clear view-field during ureteroscopic surgery, the reduction in flow may motivate an increase in inlet pressure to maintain flow and subsequent visualization. From (3.18), the inlet pressure required to maintain the steady-state flow rate at |$Q^{\textrm{SS}}$| is |$P_{\textrm{in}} = Q^{\textrm{SS}}\mu (R_{\textrm{in}}L_a+R_{\textrm{out}}L_b)$|⁠. The corresponding steady-state kidney pressure, (3.17), will simply be $$\begin{equation} P_{\textrm{k}}^{\textrm{SS}}=Q^{\textrm{SS}}\mu L_bR_{\textrm{out}}, \end{equation}$$(4.2) Fig. 6. Open in new tabDownload slide Plot (a) plots |$P_{\textrm{k}}^{\textrm{SS}}$| for a fully offset working tool as a function of working tool radius. The value of |$R_{\textrm{out}}$| is taken to be for a fully offset sheath. Plot (b) shows the effect of |$R_{\textrm{out}}$| on |$P_{\textrm{k}}$| as a function of time. Geometries corresponding to the different values are indicated. Fig. 6. Open in new tabDownload slide Plot (a) plots |$P_{\textrm{k}}^{\textrm{SS}}$| for a fully offset working tool as a function of working tool radius. The value of |$R_{\textrm{out}}$| is taken to be for a fully offset sheath. Plot (b) shows the effect of |$R_{\textrm{out}}$| on |$P_{\textrm{k}}$| as a function of time. Geometries corresponding to the different values are indicated. which is independent of |$R_{\textrm{in}}$|⁠. This indicates that if a constant flow rate is required, the kidney pressure will be independent of the size of the working tools, or their position in the working channel. Equation (4.2) is plotted as the solid black line in Fig. 6(a). In Fig. 6(b), we consider the effect of scope shaft placement within the access sheath on the time-dependent kidney pressure for fixed |$P_{\textrm{in}} = 60$| cm H|$_2$|O. We plot (3.16) (multiplied by |$P^{\star }$| to re-dimensionalize), shading the region between solutions for minimum and maximum |$R_{\textrm{out}}$| (indicated by the diagrams of offset and centred scopes). Here, we have taken |$R_{\textrm{in}}$| for an empty working channel (without a working tool). For efficient and safe ureteroscopy, we seek a combination of high flow rate and low kidney pressure. We have shown, through the wedge predictions in Figs 4(a) and (b) and 6(b), that these can be achieved simultaneously by decreasing |$R_{\textrm{out}}$|⁠. We have previously analysed optimizing the geometry of the annular sheath/scope shaft (or working tool/working channel) cross-sections to minimize axial resistance by considering elliptical cross-sections of fixed areas for the scope shaft (or working tool) and sheath (or working channel) (Williams et al., 2020). Here, we will restrict our attention to modifying the shapes of the sheath and scope shaft, as it is these that influence |$R_{\textrm{out}}$|⁠. Theoretically, we could also optimize the geometry for the working tool and working channel to minimize |$R_{\textrm{in}}$|⁠; although as discussed, this will induce the unwanted effect of decreasing |$Q$| for fixed |$P_{\textrm{in}}$|⁠. Using the optimization method outlined in Williams et al. (2020), with details provided in Appendix B, we determine the optimal shapes for the scope and access sheath, indicated in red in Fig. 6(b). This configuration has a dimensionless resistance ratio of |$\mathcal{R} = 9.6$|⁠, compared with a value of |$\mathcal{R} = 4.1$| for the offset circular scope in a circular sheath. We return to Figs 4(a) and 6(b) to show the theoretical predictions for this optimal design (dashed red lines), showing a further increase in flow rate (Fig. 4(a)), and reduction in kidney pressure (Fig. 6(b)). 5. Discussion We have considered a lumped-parameter model for irrigation during ureteroscopy. The basic modelling framework approximates the system by two distinct pipes, representing flow into the kidney through the working channel and flow back out through the access sheath, connected by a pressure point, representing the renal pressure. A similar model was previously developed (Oratis et al., 2018); although we include a data-motivated exponential, rather than linear, kidney compliance. We also improve on the model by Oratis et al. (2018), by focussing on the effects of flow resistance; in particular, modifying inflow resistance to include working tools and outflow resistance to account for the unknown position of the scope shaft within the access sheath. The chosen exponential function for kidney compliance was shown to predict a reasonable rise time of |$\mathcal{O}(10^{1}\textrm{ s})$| to steady-state flow and kidney pressure for fixed inlet pressure. This was in contrast to the predictions of a compliant linear renal stiffness (used in the model by Oratis et al., 2018), which estimated uncharacteristically long rise times of |$\mathcal{O}(10^2\textrm{ s})$|⁠. Accurately predicting the rise time to maximum renal pressure for a fixed inlet pressure is crucial. For instance, an overestimation of the rise time could lead to significant underestimates of renal pressure, leading to improper procedural decisions, such as leaving the scope inserted when withdrawal could lower pressures back to a physically safe threshold (Oratis et al., 2018). Although fitting a stiffer linear relationship between our pressure and volume also results in a predicted rise time of |$\mathcal{O}(10^1\textrm{s})$|⁠, the pressure-volume data clearly display a non-linear trend and established stiffening effects cannot be captured with a linear correlation. It is also worth noting that while an exponential function provides a good fit for pressure-volume relationship, a more detailed description would be required to make accurate mechanical predictions, for instance of the stress within the kidney tissue. Incorporating a degree of uncertainty into the model from the variable position of the scope shaft within the sheath effected predictive ‘wedges’ for steady-state flow rate and kidney pressure as functions of inlet pressure. The results of ex vivo experiments for flow rate were within the estimated range, but the renal pressure was not so accurately predicted, with measurements from one particular kidney, kidney A, falling systematically below the estimated interval. We suggested a plausible explanation for this discrepancy: variations in the position of the pressure probe within the renal cavity during the experiments. A downfall of the lumped-parameter model is that it assumes a constant renal pressure, but through numerical simulations of the steady Navier–Stokes equations in a simplified geometry, we determined that complex flow structures may exist within the kidney, leading to large spatial variations in kidney pressure. The simulations predicted a maximal variation in pressure that increased with inlet pressure and that matched, in order of magnitude at least, with the observed experimental discrepancy. While this is not proof that a different location of the pressure probe in kidney A was the cause of the discrepancy, the simulations were at least consistent with this hypothesis. The large variation in pressure in the numerical simulations also suggests that some caution may be needed when interpreting pressure readouts from the end of the ureteroscope. The presence of general terms for viscous resistance through the working channel, and back through the access sheath, provides a simple platform for theorizing on the effects of working tools and changing the geometry of the access sheath and scope shaft. We found that adding a working tool to the working channel decreases both the steady-state flow rate and kidney pressure for fixed inlet pressure. However, if a particular flow rate is sought during ureteroscopy, and the inlet pressure is adjusted accordingly to compensate for the presence of working tools, then the kidney pressure is unaffected by working tools and their size. Decreasing the outflow resistance increases flow rate and decreases kidney pressure for fixed inlet pressure. As this is a desirable combination, we determined the optimal elliptical shapes for the scope and access sheath (of fixed areas) to minimize outflow resistance. The application of mathematical models to improve ureteroscopy is emerging, with Williams (2019), Williams et al. (2019a), Williams et al. (20120), Oratis et al. (2018) and this article providing groundwork for predicting system variables and optimizing the design of pertinent biomedical devices. To bring this work closer to clinical translation, several model extensions would be useful. As we have shown, flow within the kidney can be complex, with both variable pressure and large vortical structures. These features may have a strong impact on the removal of stone dust, as well as the predictions of a systems-level model. Hence, there is a clear need to both explore flow within the kidney cavity in more detail, as well as to couple such a flow description explicitly to a systems-level description of the procedure. Similarly, our model here takes the simplest possible mechanical description of the kidney as a tissue; another useful avenue of future study would be to develop a more detailed biomechanical description. Here, a natural addition would be viscoelastic effects, as suggested by the hysteresis observed in Fig. A7. Though added modelling detail comes with a computational cost, so a significant challenge will be integrating such details within a tractable systems model. Given such advances, we can speculate on the potential integration of mathematical modelling and computational techniques in a modern operating room: we envision a ureteroscopic procedure in which the design of the medical equipment has been advised by optimization techniques; an intelligent feedback system continually monitors quantities (such as irrigation flow rates), computes hidden quantities (such as intrarenal pressure) and automatically adjusts tuneable parameters (such as inlet pressure); additionally, if a different flow rate is required, e.g. to improve visualization, predictions are made as to which factors should be varied to produce the desired result. While this describes an idealized scenario, there is a strong potential for improved procedural efficiency and reducing patient risk, and studies based on fundamental physics and mathematical modelling can provide critical first steps. Funding Engineering and Physical Sciences Research Council Centre For Doctoral Training in Industrially Focused Mathematical Modelling in collaboration with Boston Scientific (EP/L015803/1); Royal Society Leverhulme Trust Senior Research Fellowship (to S.L.W.); National Institute for Health Research Oxford Biomedical Research Centre (to B.T.). Acknowledgements The authors would like to acknowledge the support of Boston Scientific Corporation and helpful discussions with Timothy Harrah, Robert Lund, Niraj Rauniyar and Aditi Ray on the application of this work to ureteroscopy device design. The authors are also grateful for input from Simon Tavener, Patrick Farrell and Florian Wechsung on the cavity flow simulations in Section 4.2. Footnotes 1 Due to deflation of the pressure cuff, the inlet pressure slowly decreased at each ‘step’ (Fig. A7). However, the time scale was sufficiently slow that the inlet pressure could be well approximated as constant over each 5-second interval. References Amestoy , P. , Duff , I., L’Excellent , J. & Koster , J. ( 2001 ) A fully asynchronous multifrontal solver using distributed dynamic scheduling . SIAM J. Matrix Anal. Appl. , 23 , 15 – 41 . Google Scholar Crossref Search ADS WorldCat Amestoy , P. , Guermouche , A., L’Excellent , J. & Pralet , S. ( 2006 ) Hybrid scheduling for the parallel solution of linear systems . Parallel Comput. , 32 , 136 – 156 . Google Scholar Crossref Search ADS WorldCat Balay , S. , Abhyankar , S., Adams , M., Brown , J., Brune , P., Buschelman , K., Dalcin , L., Eijkhout , V., Gropp , W., Kaushik , D., Knepley , M., May , D., McInnes , L., Mills , R., Munson , T., Rupp , K., Sanan , P., Smith , B., Zampini , S. & Zhang , H. ( 2018 ) PETSc users manual . Technical Report ANL-95/11—Revision 3.9 . UChicago Argonne, LLC : Argonne National Laboratory . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Balay , S. , Gropp , W., McInnes , L. & Smith , B. ( 1997 ) Efficient management of parallelism in object oriented numerical software libraries . Proceedings of Modern Software Tools in Scientific Computing (E. Arge, A. M. Bruaset & H. P. Langtangen eds.) . Boston, MA : Birkhäuser Press , pp. 163 – 202 . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Chagnon , G. , Rebouah , M. & Favier , D. ( 2018 ) Hyperelastic energy densities for soft biological tissues: a review . J. Elasticity , 120 , 129 – 160 . Google Scholar Crossref Search ADS WorldCat Cruces , P. , Salas , C., Lillo , P., Salomon , T., Lillo , F. & Hurtado , D. ( 2014 ) The renal compartment: a hydraulic view . Intensive Care Med. Exp. , 2 , 26 . Google Scholar Crossref Search ADS PubMed WorldCat Dalcin , L. , Paz , R., Kler , P. & Cosimo , A. ( 2011 ) Parallel distributed computing using Python . Adv. Water Resour. , 34 , 1124 – 1139 . Google Scholar Crossref Search ADS WorldCat Deuflhard , P. ( 2011 ) Newton Methods for Nonlinear Problems . Berlin, Germany : Springer . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Farrell , P. , Mitchell , L. & Wechsung , F. ( 2019 ) An augmented Lagrangian preconditioner for the 3D stationary incompressible Navier–Stokes equations at high Reynolds number . SIAM J. Sci. Comput. , 41 , A3073 – A3096 . Google Scholar Crossref Search ADS WorldCat Gregersen , H. ( 1996 ) Biomechanics of the Gastrointestinal Tract , vol. 4 . London : Springer . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Griffiths , D. J. ( 1980 ) Urodynamics: The Mechanics and Hydrodynamics of the Lower Urinary Tract . Bristol : Hilger . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hendrickson , B. & Leland , R. ( 1995 ) A multilevel algorithm for partitioning graphs . Supercomputing ’95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM) . New York : ACM Press , p. 28 . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Heyda , J. ( 1959 ) A Green’s function solution for the case of laminar incompressible flow between non-concentric circular cylinders . J. Franklin Inst. B , 267 , 25 – 34 . Google Scholar Crossref Search ADS WorldCat Holzapfel , G. ( 2001 ) Biomechanics of soft tissue . Handbook of Material Behavior Models . Volume III, Multiphysics Behaviors, Chapter 10. (Composite Media ed.). Boston : Academic Press , pp. 1049 – 1063 . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Humphrey , J. ( 2003 ) Continuum biomechanics of soft biological tissues . Proc. Royal Soc. A , 459 , 3 – 46 . Google Scholar Crossref Search ADS WorldCat Kauer , M. ( 2001 ) Inverse finite element characterization of soft tissues with aspiration experiments . Ph.D. Thesis , ETH Zürich , Swiss Federal Institute of Technology Zürich. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kirby , R. & Mitchell , L. ( 2018 ) Solver composition across the PDE/linear algebra barrier . SIAM J. Sci. Comput. , 40 , C76 – C98 . Google Scholar Crossref Search ADS WorldCat Kum , F. , Mahmalji , W., Hale , J., Thomas , K., Bultitude , M. & Glass , J. ( 2016 ) Do stones still kill? An analysis of death from stone disease 1999–2013 in England and Wales . BJU Int. , 118 , 140 – 144 . Google Scholar Crossref Search ADS PubMed WorldCat Lamb , H. ( 1895 ) Hydrodynamics . Oxford : Jesus College, Cambridge University Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Mitchell , L. & Müller , E. ( 2016 ) High level implementation of geometric multigrid solvers for finite element problems: applications in atmospheric modelling . J. Comput. Phys. , 327 , 1 – 18 . Google Scholar Crossref Search ADS WorldCat Oratis , A. , Subasic , J., Hernandez , N., Bird , J. & Eisner , B. ( 2018 ) A simple fluid dynamic model of renal pelvis pressures during ureteroscopic kidney stone treatment . PLoS One , 13 , e0208209 . Google Scholar Crossref Search ADS PubMed WorldCat Rathgeber , F. , Ham , D., Mitchell , L., Lange , M., Luporini , F., McRae , A., Bercea , G.-T., Markall , G. & Kelly , P. ( 2016 ) Firedrake: automating the finite element method by composing abstractions . ACM Trans. Math. Softw. , 43 , 24:1 – 24:27 . OpenURL Placeholder Text WorldCat Shattock , S. ( 1905 ) A Prehistoric or Predynastic Egyptian Calculus . Egypt : Adlard . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wechsung , F. ( 2019 ) Shape optimisation and robust solvers for incompressible flow (Chapter 7) . Ph.D. Thesis , University of Oxford , Bodleian Libraries, Oxford . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Williams , J. ( 2019 ) Mathematical modelling of fluid flows during ureteroscopic kidney stone removal . Ph.D. Thesis , University of Oxford , Bodleian Libraries, Oxford . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Williams , J. , Turney , B., Moulton , D. & Waters , S. ( 2020 ) Effects of geometry on resistance in elliptical pipe flows . J. Fluid Mech. , 891 , A4 . Google Scholar Crossref Search ADS WorldCat Williams , J. , Turney , B., Rauniyar , N., Harrah , T., Waters , S. & Moulton , D. ( 2019a ) The fluid mechanics of ureteroscope irrigation . J. Endourol. , 33 , 28 – 34 . Google Scholar Crossref Search ADS WorldCat Wilson , W. T. & Preminger , G. M. ( 1990 ) Intrarenal pressures generated during flexible deflectable ureterorenoscopy . J. Endourol. , 4 , 135 – 141 . Google Scholar Crossref Search ADS WorldCat Appendix. A. Data analysis A.1 Raw data The independent variable in the experiments, |$P_{\textrm{in}}$|⁠, is shown as a function of time for each experimental run in Fig. A7. The three rows show data for kidneys A, B and C, respectively, and the corresponding data set (I or II) is indicated in the figure titles. A.2 Data set I Intrarenal pressure is plotted against mass for data set I, (Fig. A7(e), (h), (i) and (j)) in Fig. A8(a). The mass and pressure data have been appropriately interpolated so they can be plotted against each other, with time parameterizing the curve. The interval including both increasing and decreasing |$P_{\textrm{in}}$| is plotted, providing both ‘loading’ and ‘unloading’ curves. We note the presence of hysteresis, particularly prevalent for the data from kidney A in Fig. A8(a), i.e. the loading and unloading data follow different curves. Additionally, with repeated loading and unloading cycles on a single kidney (kidney C), the hysteretic curves shift to the right. Fig. A7. Open in new tabDownload slide Inlet pressure (⁠|$P_{\textrm{in}}$|⁠) over time. All vertical axes go from 0-500 cm H|$_2$|O and horizontal axes from 0-400 s. Fig. A7. Open in new tabDownload slide Inlet pressure (⁠|$P_{\textrm{in}}$|⁠) over time. All vertical axes go from 0-500 cm H|$_2$|O and horizontal axes from 0-400 s. Fig. A8. Open in new tabDownload slide All continuous-loading data sets. |$P_{\textrm{k}}$| is plotted against |$M_{\textrm{k}}$| for the loading and unloading curves for (e), (h), (i) and (j) in Fig. A7. Fig. A8. Open in new tabDownload slide All continuous-loading data sets. |$P_{\textrm{k}}$| is plotted against |$M_{\textrm{k}}$| for the loading and unloading curves for (e), (h), (i) and (j) in Fig. A7. Fig. A9. Open in new tabDownload slide All continuous-loading data sets. |$P_{\textrm{k}}$| is plotted against |$M_{\textrm{k}}$| for the loading and unloading curves for (e), (h), (i) and (j) in Fig. A7. Fig. A9. Open in new tabDownload slide All continuous-loading data sets. |$P_{\textrm{k}}$| is plotted against |$M_{\textrm{k}}$| for the loading and unloading curves for (e), (h), (i) and (j) in Fig. A7. It is well established that biological tissue displays properties of viscoelastic solids (Humphrey, 2003). Thus, stress depends not only on the applied strain but also on the rate of strain. This leads to viscoelastic features such as hysteresis (difference in loading and unloading stress-strain curves), and stress relaxation, where if a material is suddenly strained and the strain is maintained constant, the stress in the material will decrease over time (Gregersen, 1996). Additionally, with repeated loading and unloading cycles, the load-deformation curves will shift to the right, and the hysteretic effects will diminish. By repeated cycling, eventually a steady-state is reached at which no further change will occur unless the cycling routine is changed. In this state the tissue is said to be preconditioned (Kauer, 2001). By looking at the plots of inlet pressure over time (Fig. A7), we see that Fig. (h–j) will lend themselves best to tissue preconditioning. This is reflected in Fig. A8(a), where we see that the corresponding curves shift to the right, and curves (i) and (j) appear to be nearly indistinguishable. We thus use the loading curves of these two data sets in Section 4.1 to estimate the kidney compliance. We note, however, that all loading curves fit the exponential trend of (3.3) shown as the red lines in Fig. A8(b). This further validates the use of an exponential kidney compliance, where |$P^{\star }$| and |$C$| may depend on the kidney, loading rate and any previous loading. We have chosen two data sets for kidney C to determine the best-fit values for the exponential compliance model |$P^{\star } = 5.24\textrm{ cm H}_2\textrm{O}$| and |$C = 0.43\textrm{ mL}^{-1}$|—(4.1)—as for these runs we could be most confident that the kidney had been preconditioned. We note, however, that choosing the loading curve for kidney A (the black symbols in Fig. A8) results in values of |$P^{\star } = 8.83\textrm{ cm H}_2\textrm{O}$| and |$C = 0.35\textrm{ mL}^{-1}$|⁠. The data for kidney A (black triangles) along with the best-fit exponential curve (red line) are shown in Fig. A9. Using these values in the model for kidney pressure, (3.16), we obtain the red curve in Fig. A9(b). As a comparison, we reproduce the red curve from Fig. 3(b), which uses the fitted |$P^{\star }$|⁠, |$C$| values for kidney C as the dashed black line in Fig. A9(b). Of interest is to note the similarity of the two curves produced with different values of |$P^{\star }$| and |$C$|⁠, indicating that different kidney compliances will yield similar qualitative—and even quantitative—kidney pressures over time. A.3 Data set II The data in set II, with inlet pressure shown in Fig. A7(a–d), (f), (g), (j) and (k), is used to estimate steady-state flow rate and pressure data. Multiple 5-second intervals, where the inlet pressure is approximately constant, are extracted from each data set. The flow rate as a function of time is estimated by taking the difference between successive outflow mass measurements, divided by the time step. A plot of the mean flow rate against the mean inlet pressure over each 5-second interval provides the data points in Fig. 4(a). The error bars in Fig. 4(a and b) combine the standard deviation between successive mass measurements within each 5-second interval and the measurement error, as described in the following section. Note that all data from set II is included in Fig. 4; since the data from different runs on each kidney were very similar, we have not delineated between the runs, only between the different kidneys. A.4 Error estimation The mass balance used was the OHAUS Traveler, which has a readability of 0.01 g. Therefore, the systemic error in mass measurements was taken to be |$\varDelta _M = \pm 0.005$| g. This is displayed as the (indiscernible) horizontal error bars in Figs 3, A8 and A9. A FISO Fiber Optic Pressure Sensor FOP-M260-21 was used to record pressure measurements and has an accuracy of |$0.5\%$| of the |$-300$| to |$300$| mmHg pressure scale. Therefore, the systematic error in pressure measurements was taken to be |$\varDelta _P = \pm 0.005\times 600 = \pm 3$| mmHg. This is taken as the vertical error bars in Figs 3, A8 and A9 (conversion between mmHg and H|$_2$|O requires multiplication by a factor of approximately |$1.36$|⁠). For the steady-state experiments, we calculate a series of |$n$| flow rates over a 5-second interval $$\begin{equation} Q_i = {\sum_{i=1}^{n}}\frac{M_{\textrm{k}}^{i+1}-M_{\textrm{k}}^{i}}{\varDelta t_i}, \end{equation}$$(A.1) where |$\varDelta t_i = 500$| ms (so |$n \approx 10$|⁠). The standard deviation of the series of flow rates then has an associated mean |$\mu _Q$| and standard deviation |$\varsigma _Q$|⁠. As flow rate is estimated by subtracting successive mass measurements, the total systematic error in flow rate was |$\varDelta _{Q} = \sqrt{2}\varDelta _M\sqrt{\sum _{i=1}^n(1/\varDelta t_i)^2}/(n-1)$|⁠. We assume that the standard deviation in the flow rates and the measurement errors due to the readability of the mass balance are independent and thus the total error $$\begin{equation} T.E. = \sqrt{\varsigma_Q^2 + \varDelta_Q^2}, \end{equation}$$(A.2) is taken as the vertical error bars in Fig. 4(a). The vertical error bars on the approximated steady-state pressure in Fig. 4(b) is given by $$\begin{equation} T.E. = \sqrt{\varsigma_P^2 + \varDelta_P^2}, \end{equation}$$(A.3) where |$\varsigma _P$| the standard deviation of the approximately |$625$| pressure measurements over each 5-second interval. Appendix. B. Determining the optimal elliptical shape In our previous study of flow through elliptical annular pipes (Williams et al., 2019a), we sought to determine the optimal elliptical geometry—elliptical shapes of the outer and inner boundaries and the position and orientation of the annulus core—to minimize flow resistance. We formulate the optimization problem $$\begin{equation} \min_{\boldsymbol{g}} -Q, \qquad \text{s.t.}\ \boldsymbol{c} \geqslant 0, \end{equation}$$(B.1) where |$Q$| is determined as the solution to $$\begin{equation} \nabla^2w = -1, \qquad Q = \iint_{\Omega}w\:\textrm{d}\Omega, \end{equation}$$(B.2) and |$\varOmega $| is the annular space between the inner and outer ellipses representing the scope shaft and access sheath. In (B.1), |$\boldsymbol{c}$| is a vector of six geometric non-linear constraints which ensure that the inner ellipse is fully contained within the outer ellipse (Williams et al., 2019a). We solve equation (B.2) using open-source finite element library Firedrake and implement MATLAB’s optimization routine, fmincon, to solve the optimization problem (B1) with an interior-point method from eight randomly generated starting points in the feasible 5D parameter space. As the optimization routine seeks local minima, a multi-start approach allows for global optimization. The optimal geometry found through this method is shown in Fig. B10 with major and minor axis of the scope shaft |$(0.17, 0.15)$| cm and of the access sheath |$(0.23, 0.16)$| cm indicated. The contour lines in Fig. B10 provide contours of |$w$|⁠; |$Q \approx 1.2$| is obtained by (B2)b. Fig. B10. Open in new tabDownload slide The optimal geometry for scope shaft/access sheath with measurements provided in cm. The colourbar provides contours of |$w$| (dimensionless). Fig. B10. Open in new tabDownload slide The optimal geometry for scope shaft/access sheath with measurements provided in cm. The colourbar provides contours of |$w$| (dimensionless). Appendix. C. Pressure variation in the renal pelvis To model pressure variation in the renal pelvis (Fig. 5), we consider the light-pink 2D domain in Fig. C11, where irrigation fluid enters through |$\varGamma _{\textrm{in}}$| (the tip of the working channel) and exits through |$\varGamma _{\textrm{out}}$| (the gap between the scope shaft and the access sheath. Fig. C11. Open in new tabDownload slide The geometric domain where |$a$| represents the radius of the working channel, |$b$| the radius of the scope shaft minus the radius of the working channel and |$c$| the gap between the scope and access sheath for a centred scope. The light pink domain represents the domain, |$\varOmega _{\textrm{k}}$|⁠, in which the Navier–Stokes equations are solved numerically. Fig. C11. Open in new tabDownload slide The geometric domain where |$a$| represents the radius of the working channel, |$b$| the radius of the scope shaft minus the radius of the working channel and |$c$| the gap between the scope and access sheath for a centred scope. The light pink domain represents the domain, |$\varOmega _{\textrm{k}}$|⁠, in which the Navier–Stokes equations are solved numerically. We model irrigation flow in the cavity with the steady Navier–Stokes equations, considering a velocity field in Cartesian coordinates, |$\mathbf{u} = u\mathbf{i} + v\mathbf{j}$|⁠. We impose a fully developed parabolic profile with maximum velocity |$U$| at the inlet boundary. Non-dimensionalizing |$(x, y)$| with the half-width of the inflow channel (⁠|$a$| in Fig. C11), |$\mathbf{u}$| with |$U$|⁠, and pressure |$p$| with |$\mu U/a$|⁠, the dimensionless steady Navier–Stokes equations are $$\begin{equation} \textrm{Re}(\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p +\nabla^2\mathbf{u}, \qquad \nabla\cdot\mathbf{u} = 0, \end{equation}$$(C.1) where |$\textrm{Re} = \rho Ua/\mu $|⁠. The prescribed dimensionless inflow is then $$\begin{equation} u = 1-y^2, \qquad v = 0, \qquad \textrm{on } \varGamma_{\textrm{in}}. \end{equation}$$(C.2) At the outflow boundaries |$\varGamma _{\textrm{out}}$|⁠, we prescribe zero normal and tangential stress, |$\boldsymbol{\sigma }\mathbf{n} = \boldsymbol{0}$|⁠, which provides $$\begin{equation} 2u_x - p =0, \qquad u_y + v_x = 0, \qquad \textrm{on }\Gamma_\textrm{out}. \end{equation}$$(C.3) On the walls (black solid lines in Fig. C11), we assume no-slip $$\begin{equation} u = v = 0, \qquad \textrm{on}\ \varGamma_{\textrm{wall}}. \end{equation}$$(C.4) For a particular experiment of flow rate |$Q$|⁠, we can approximate the maximum inflow velocity, |$U = 2(Q/\pi a^2)$|⁠, and hence the Reynolds number $$\begin{equation} \textrm{Re} = \rho2Q/\pi a\mu. \end{equation}$$(C.5) Equations (C1) are solved subject to (C2) and (C3), and (C4) are solved with finite elements implemented in Firedrake (Amestoy et al., 2001, 2006; Balay et al., 2018, 1997; Dalcin et al., 2011; Deuflhard, 2011; Hendrickson & Leland, 1995; Kirby & Mitchell, 2018; Mitchell & Müller, 2016; Rathgeber et al., 2016). The solver we use was provided by the authors of Farrell et al. (2019) but extends the work of that paper to use the exactly divergence-free Scott–Vogelius element instead, as explained in (Wechsung, 2019). The relevant Reynolds numbers for simulations representing the flow rates in Fig. 4(a), and thus measured renal pressures in 4(b), range from |$\textrm{Re} = 0$| to |$\textrm{Re} = 1883$|⁠. Results of the simulation for |$\textrm{Re} = 1500$| (⁠|$Q \approx 1.4$| cm|$^3$|/s) are shown in Fig. 5(a and b). We have previously determined the existence of multiple solutions to the steady Navier–Stokes equations solved in the domain pictured in Fig. C11, with the asymmetric solution branch theorized to be stable, where it exists (Williams, 2019a). By initially placing a no-slip condition on one of the outflow boundaries, we ensure the asymmetric solution branch. © The Author(s) 2020. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A lumped-parameter model for kidney pressure during stone removal
JF - IMA Journal of Applied Mathematics
DO - 10.1093/imamat/hxaa020
DA - 2020-09-25
UR - https://www.deepdyve.com/lp/oxford-university-press/a-lumped-parameter-model-for-kidney-pressure-during-stone-removal-WvoZyaJyUx
SP - 703
EP - 723
VL - 85
IS - 5
DP - DeepDyve
ER -