Biomechanical modelling of bile flow in the biliary system

The biliary system consists of the biliary tree, gallbladder and major duodenal papilla. Soft tissues compliance plays important role in the bio-fluids transport. Particularly, bile flow disturbances due to bile duct wall motor function changes in the extra-hepatic ducts, from medicine point of view are called dyscinesia of biliary tract. Fluid motion in the elastic and compliant ducts can be described by different models (for example, Windkessel model, peristaltic fluid motion, FSI algorithm). Our approach is decomposition of the biliary system into three compartments (extra-hepatic biliary tree, gallbladder, major duodenal papilla). Bile flow in the extra-hepatic ducts is simulated using FSI algorithm. Bile flow in the gallbladder can be described as flow in the reservoir with compliant ducts using Windkessel model. Bile flow in the major duodenal papilla is considered as peristaltic fluid motion, because the wall contraction is really important factor of fluid motion in that segment. The coupling of these compartments is performed by boundary conditions. The biliary system geometry was obtained using MRI patient-specific data. It was confirmed that normal bile can be modeled as Newtonian fluid and lithogenic bile can be modeled as non-Newtonian fluid (Carreau fluid). Bile ducts were modeled as hyperelastic material.


Introduction
The biliary system is responsible for bile producing, storing and transport to the duodenum in order to emulgate fats.Despite several biomechanical models of bile flow in different segments, there is no entire biomechanical model.Our approach is related to considering of the biliary system as three compartments.Each of compartment is modelled using different models depending on mechanical performance of this compartment based on physiological data.Computational models of bile flow dynamics were proposed by Luo et al [1], [2], [3], Agrawal et al [4], [5], Kuchumov et al [6], [7], Al-Atabi et al [8], Ro-Chin Lo et al [9].Peristaltic bile flow models were proposed by Maiti and Misra [10] and Kuchumov et al [11], [12].There is no entire biomechanical model of the bile flow in the biliary system and the purpose of this paper is present a such model.
Our approach is decomposition of the biliary system into three compartments (extrahepatic biliary tree, gallbladder, major duodenal papilla) using several models and approaches.The compartments coupling into the entire model is performed using boundary conditions.
• Bile flow in the gallbladder.
We used a modified Windkessel model to evaluate bile flow in the patient-specific gallbladder model.Using ultrasound data we can evaluate Windkessel model parameters for the patient and non-invasively determine pressure changes on time during the gallbladder refilling and emptying.The pressure on time dependence is applied as a boundary condition when we investigate bile flow in the biliary tree.
• Bile flow in the biliary tree.
We considered 4 cases: healthy bile flow (considered as Newtonian fluid), lithogenic bile (considered as non-Newtonian fluid (Carreau fluid)) flow, lithogenic bile flow in the case of gallstone presence; bile flow after cholecystectomy (the gallbladder removal).The velocity and pressure distributions during the gallbladder emptying were presented.Velocities were found to have lower magnitude in the case of lithogenic bile, but the pressures are higher in this case.Stress state of bile ducts was also evaluated.It was shown that von Mises stress distribution is located in the cystic duct and hepatic ducts mostly; whereas, the shear stress distribution is mostly prevailing in the common hepatic duct.It is believed that duct wall shear stress can be considered as indicator of gallstones formation and von Mises stress is related to biliary pain.The developed model can be applied in the medical practice to evaluate the circumstances of surgical interventions.
• Bile flow in the major duodenal papilla.
Choledochopancreatic reflux (i.е. the flow of the gallbladder bile coming out the common bile duct into the pancreatic duct instead of the duodenum) is known to be one of the reasons of the pancreatitis (inflammation of the pancreas).Understanding of the reasons of the reflux from physiological, hydrodynamic, biomechanical points of view is still a challenging task.The current paper aims at developing mathematical model of the peristaltic bile transport flow through the duct at papillary stenosis as a tapered finite-length tube.It allows evaluating velocities and pressure distribution along the tube, and detecting choledochopancreatic reflux occurrence conditions.Adopting the perturbation method, the analytical solutions for velocities and pressures are obtained.Pressure distribution versus axial coordinate at different time instants are plotted for various values of Weissenberg number and amplitude ratio.It revealed that the amplitude ratio has more effect on the pressure distribution along the tube compared to the Weissenberg number.The values of the pressure gradient corresponding to reflux occurring are obtained.Moreover, it is reported that the pressure drop value corresponding to average flow rate equal to zero may serve as reflux occurrence criterion.during emptying scheme .

Modified Windkessel model
The flow scheme due to pressure gradient and interaction with elastic walls can be described by Windkessel model [14]: where V is the gallbladder volume, С is a gallbladder wall compliance, Q is a bile flow rate, p is a pressure inside the gallbladder, pd is a pressure in the duodenum (pd = 6 mm Hg) [15], R is a hydraulic resistance.Substituting (3) in (1) with account of (2), considering С = const, one can obtain: where pe is the pressure when gallbladder is empty (pe = 11 mm Hg) [12], te is the emptying time.Volume change dependence :

Ultrasound gallbladder performance
In order to find patient-specific parameters of Windkessel model, the experiment was performed.The ultrasound measurements of volunteer's full gallbladder volume were made initially.The typical ultrasound image of the gallbladder is presented in Fig. 4.After that, the volunteer had a meal in order to force the gallbladder emptying.After 20 minutes when the meal was taken with 5 minute interval, the volume measurements were made.The measurements were made from 20 to 240 minutes after taking a meal.After that using least square method, we obtained Windkessel model parameters can be found by (6) (Fig. 5).Using the obtained parameters, we can substitute them into the pressure on time dependence and use it as a boundary condition, when bile flow modelling in the biliary tree is performing.

Bile flow in the biliary tree
The second stage of simulation of bile flow in the biliary system is FSI modeling of the bile flow in the biliary tree.To make a patient-specific FSI modeling you need to create 3D patient-specific geometry; to evaluate fluid and solid properties Symmetric finite element model was created using ANSYS APDL Mechanical (ANSYS, Inc., Canonsburg, PA, USA).The dimensions of the model were taken from experiment images processing using ImageTool software (UTHSCSA, USA).The thickness was taken as 10 -3 m.
The obtained parameters were exported to the ANSYS Workbench 13.0 (ANSYS, Inc., Canonsburg, PA, USA) (Engineering Data), which allowed to recompute stress-strain curves for uniaxial, biaxial and shear experiments for the hyperelastic material described by Mooney-Rivlin 2 parameter model.

Fluid-Structure Interaction Problem Statement
Mass and momentum conservation equations for an incompressible fluid can be expressed as where is f  is the fluid density, p is the pressure, u is the fluid velocity vector and ug is the moving coordinate velocity.In the Arbitrary Lagrangian -Eulerian formulation, (u−ug) is the relative velocity of the fluid with respect to the moving coordinate velocity.Here τ is the deviatoric stress tensor.This tensor is related to the strain rate tensor; if the fluid is incompressible and viscosity is constant across the fluid, this equation can be written in terms of an arbitrary coordinate system as where xj is the th spatial coordinate, vi is the fluid's velocity in the direction of axis i, τij is the jth component of the stress acting on the faces of the fluid element perpendicular to axis.
For healthy bile, which is considered as Newtonian Fluid const   , The Carreau's equation, which is used for simulation of lithogenic bile, is written as The momentum conservation equation for solid body can be written as where s  , s  , and g u are solid density, solid tensor, and local acceleration of the solid, respectively.
It is known that blood vessels can be described as hyperelastic materials [16].Due to similar anatomical composition, the bile ducts can also be considered as hyperelastic materials.For hyperelastic materials, the stress-strain relation is written as follows The Mooney-Rivlin hyperelastic potential takes on the form The FSI interface should satisfy the following conditions: 1.The displacements of the fluid and solid domain should be compatible i.e.
2. The tractions at this boundary must be at equilibrium: 3. and the no-slip condition for the fluid should satisfy In the above conditions,  ,  and n are respectively, displacement, stress tensor and boundary normal with the subscripts 'f' and 's' indicating a property of the fluid and solid.

Solid Domain
As shown in Fig. 6, the ends of the extra-hepatic biliary tree were constrained by specifying zero-displacements in all directions and prohibition of the rotation about all axes (boundary condition type: fixed support).

Fluid domain
During the gallbladder refilling, the bile comes out from the the liver via right and left hepatic ducts and flows into the gallbladder.Bile pressures and initial velocity were applied to the inlets and the outlet (Fig. 6).Following the work of Howard et al [17], who measured the mean flow rate of the bile from the gallbladder (2.0-3.0 ml/min) after meal.Using relation where v is velocity, S is the cross section area, which can be evaluated for the given patientspecific model, we expressed the velocity from this equation and obtained the mean value of velocity 6 mm/s and used this value as inlet velocity.Outlet pressure corresponds to p(t) obtained in the previous section.Pressure in the end of the biliary tree was taken equal to the pressure in the duodenum i.e. 0.96 kPa [18].Fluid and solid meshes consisted of approximately 610,000 and 79,000 elements, respectively.Solid mesh was mapped with tetrahedral element size 0.3 mm.The fluid mesh was free with tetrahedral elements also.One-way FSI algorithm was adopted and the simulations were performed using ANSYS Workbench 13.0 (ANSYS, Inc., Canonsburg, PA, USA).The two solvers (ANSYS Mechanical and ANSYS CFX) were coupled and solved iteratively.The time step used for this simulation was 0.04 s.Initially, the fluid field is solved until the convergence criteria are reached.The calculated forces at the structure boundaries are then transferred to the structure side.The structure side is calculated until the convergence criterion is reached.The solution is finished when the maximum number of time steps (100) is reached.

Velocity distributions
Velocity distributions for two cases (flow of lithogenic bile in the biliary tree and bile flow in the extra-hepatic tree after cholecystectomy) are presented in Fig. 7.When gallbladder empties, the bile flows out of the gallbladder via cystic duct and common bile duct.At the same time bile flows out of the liver via hepatic ducts and common bile duct.We simulated these physiological phenomena in the present investigation.From Fig. 7 one can see that the peak of velocity is observed in the cystic duct, where bile accelarates towards the common bile duct.In our opinion, it is related with the fact that cystic duct has very complex inner geometry (spiral valves of Heister), which function is to damp such high velocities.In this case only external geometry was created without account of inner structure.As it was mentioned previously, the bile flow velocities play role in the stone formation.Lower velocity fields are potential zones for stone occurence.The minimal velocities are observed in the case of bile flow in the biliary tree after cholecystectomy.It is related with the fact that the gallbladder and cystic duct (where commonly the maximal velocity is observed) are absent.The bile amount is less and thus the velocity are lower than in the other cases.It also stated that lower velocities may lead to the stone formation.Our results showed that after cholecystectomy bile velocity significantly decreases and stone formation is probable for the current situation.The comfirmed fact of repeated stone formation after cholecystectomy can be also found in the medical literature [19].

Bile flow in the major duodenal papilla
We consider the peristaltic flow of a non-Newtonian fluid (Carreau's fluid) in a tapered finite-length tube.Sinusoidal waves of constant speed propagate along the channel boundaries.The wavelength is comparable with the channel length (L≈λ), thus the wave number is very small, and Reynolds number is negligible.A wall geometry is described by the relation where a is an inlet radius; l is a slope coefficient; c is a wave velocity; b is an amplitude of peristaltic wave; Z is a longitudinal coordinate.
The equations of motion are where , W U are axial and radial velocities.Let us introduce non-dimensional variables: After non-dimensionalization, the equations are as follows   To reduce the problem, we will use the lubrication theory approximations of infinitesimally small wall curvature (δ→0) and small Reynolds number (Re→0).The approximations assume that inertial effects are negligible and that the dominant axial scale is much larger than the dominant radial scale [20]. Thus, where Boundary conditions are written as follows: ----for axial velocity ----for radial velocity The finite tube length requires the pressures at two ends of the tube In order to solve equations ( 24)-( 26) and account for boundary conditions ( 27)-(32), the perturbation method is adopted and the solution is represented as an expansion over Weissenberg number (We): where 0 ( ) where is equal to

Influence of Weissenberg number on the pressure distribution
This section is devoted to the analysis of Weissenberg number (We) and amplitude ratio (φ) influence on the pressure distribution [p(ξ, t) − p(0, t)] along the tube, and dependence of averaged flow rate (Q) versus pressure to analyze the reflux occurring conditions.Figures are plotted by using MatLab software.
Pressure distributions at various Weissenberg numbers (We = 0.1; 0.5) are shown in Figs.8-9.The pressure distribution evolution along the tapered tube length at Weissenberg number We = 0.5 is displayed in Figs. 9, a-f.Fig. 9, a (t = 0) shows that at first the pressure decreases slowly and finally rises sharply to meet the leading end of the bolus.After one half of the periodic cycle (Fig. 9, d), pressure value reaches its minimum.Finally, at the time t = 1 (Fig. 9, f), which represents the completion of one period, the pressure distribution resembles the unity at t = 0; it indicates that a new cycle is about to begin (Fig. 9, a).The same behavior can be seen in Fig. 5 plotted for Weissenberg number We = 0.9.It is observed that pressure difference in different segments of the tube increases with the increase of Weissenberg number.It can be physically interpreted that if the fluid viscosity is less, then less pressure is required for the Carreau's fluid propagation through the tapered tube.

Influence of amplitude ratio on the pressure distribution
Influence of amplitude ratio on the pressure distribution along the tube is shown in Figs.
10-11.The calculations were made for three amplitude ratios: φ = 0.1; 0.2.Weissenberg number value was taken equal to 0.5.The increase of amplitude ratio decreases crosssection during bolus propagation and simultaneously increases the pressure value difference.It can be seen also that amplitude ratio has more influence on the pressure distribution compared with the Weissenberg number.

Conclusions
According to [21], the use of elective cholecystectomy has increased dramatically following the widespread adoption of laparoscopic cholecystectomy (Fig. 6).The number of complications is still high.One of the reasons is the absence of patient-specific approach for patient's treatment.The surgeons use their own experience, but not use the results of biomechanics to evaluate the intervention outcomes.Our research will help surgeons to decrease the number of complications in population and will shorten the recovery time after surgery.
The reported study was funded by RFBR according to the research project № 16-08-00718 А and Russian Government Assignment grant 19.7286.2017/8.9.

Fig. 7 .
Fig. 7. Velocity distributions in case of lithogenic bile flow and bile flow in the extra-hepatic tree after cholecystectomy.
Solving the system (24)-(26) subjected to the boundary conditions, the axial and radial velocities are obtained in the following form: