Med Biol Eng Comput (2014) 52:841–849 DOI 10.1007/s11517-014-1187-1
ORIGINAL ARTICLE
Backflow length predictions during flow‑controlled infusions using a nonlinear biphasic finite element model Gustavo A. Orozco · Joshua H. Smith · José J. García
Received: 20 November 2013 / Accepted: 13 August 2014 / Published online: 26 August 2014 © International Federation for Medical and Biological Engineering 2014
Abstract A previously proposed finite element model that considers geometric and material nonlinearities and the free boundary problems that occur at the catheter tip and in the annular zone around the lateral surface of the catheter was revised and was used to fit a power-law formula to predict backflow length during infusions into brain tissue. Compared to a closed-form solution based on linear elasticity, the power-law formula for compliant materials predicted a substantial lower influence of the shear modulus and catheter radius on the backflow length, whereas the corresponding influence for stiffer materials was more consistent with the closed-form solution. The finite element model predicted decreases of the backflow length for reduction of the shear modulus for highly compliant materials (shear modulus less than 500 Pa) due to the increased area of infusion and the high fluid fraction near the infusion cavity that greatly increased the surface area available for fluid transfer and reduced the hydraulic resistance toward the tissue. These results show the importance of taking into account the material and geometrical nonlinearities that arise near the infusion surface as well as the change of hydraulic conductivity with strain for a proper G. A. Orozco Escuela de Ingeniería Mecánica, Universidad del Valle, Cali, Colombia e-mail: [email protected] J. H. Smith (*) Department of Mechanical Engineering, Lafayette College, Easton, PA 18042, USA e-mail: [email protected] J. J. García Escuela de Ingeniería Civil y Geomática, Universidad del Valle, Cali, Colombia e-mail: [email protected]
characterization of backflow length during flow-controlled infusions into the brain. Keywords Convection-enhanced delivery · CED · Infusion · Biphasic materials
1 Introduction Convection-enhanced delivery (CED) is a means to bypass the blood–brain barrier and deliver controlled volumes of drugs over localized zones of the brain for the treatment of diseases and tumors [3, 19]. Many experimental studies have been performed in agarose gels and animals to understand the effect of different parameters on the outcome of CED. Early brain infusion experiments with animals at higher flow rates showed backflow, in which an annular zone is formed outside the catheter and the infused drug preferentially flows toward the surface of the brain rather than through the tissue toward the area targeted for delivery [6]. In this case, the transport of the infused drug mainly occurs from the annular zone around the catheter rather than from its tip. To provide recommendations for clinical studies, controlled experiments using agarose gels as a brain tissue phantom have been performed for the measurement of infusion pressure, drug distribution, and backflow length [8, 9]. The foregoing experimental and clinical studies have been complemented with analytical models aimed to study the variables that affect backflow. The seminal model by Morrison et al. [18], later enhanced by Raghavan et al. [20], allows for the prediction of the backflow length and fluid pressure by considering the deformation of the tissue around the external boundary of the catheter, the axial flow in the annular gap developed around the cannula, and the
13
842
Med Biol Eng Comput (2014) 52:841–849
radial flow from this annular region into the porous tissue. Specifically, Morrison et al. [18] developed a power-law correlation for the backflow length
L ∝ Q0.6 R0.8 rc0.8 G−0.6 µ−0.2 , 0
(1)
where Q is the infusion flow rate, R is a tissue hydraulic resistance, rc is the catheter radius, G0 is the tissue shear modulus, and µ is the fluid viscosity. However, this formula was derived under some limiting assumptions, such as considering the solid phase of the infused tissue as a linearly elastic material under infinitesimal deformations, whereas mechanical testing has shown large deformations under physiological loadings [10, 15–17]. As recommendations for defining the parameters of infusions during clinical studies are being influenced by the backflow predictions of the analytical models [21–23] and as experimental infusion studies in animals have shown that deformation of the brain may be significant [28], it is important to develop models that incorporate the nonlinear behavior of brain tissue. To relax the assumptions of the analytical solution and to include geometric and material nonlinearities documented in experiments, we developed a finite element (FE) model of backflow in brain tissue phantoms [11] that was validated with published experimental results [20]. In this study, we updated our finite element model of backflow with a material model for brain tissue that is capable of reproducing the nonlinear stress–strain behavior that has been documented experimentally [10, 15–17]. We then used our revised finite element model to develop a powerlaw type formula similar to Eq. (1) that relates backflow length to the infusion flow rate and properties of the infused tissue. Lastly, we assessed the influence of nonlinear effects in the development of backflow.
2 Methods 2.1 Geometry and boundary conditions Our FE model of infusion into brain tissue from catheters of 0.105, 0.33, 0.50, and 0.98-mm radius was developed using ABAQUS 6.9-1 (Simulia, Providence, RI) considering axial symmetry and the planar geometry of Fig. 1. The computational domain extended between the outer radius of the catheter and a cylinder of 20-mm radius, with the lower portion of the model represented as a 20-mm radius hemisphere, consistent with previous studies [11, 18, 20]. Following a mesh convergence study based on the maximum fluid pressure, the domain was meshed with 4,486, 3,431, 3,364, and 3,241 elements of type CAX4P for the 0.105, 0.33, 0.50, and 0.98-mm radius catheters, respectively. These elements form an axisymmetric domain with respect to the vertical z-axis.
13
Fig. 1 Sketch of the domain considered in the FE model of backflow around a 0.98-mm radius catheter
Fig. 2 Representation of the radial displacement ur, radial fluid velocity qr, fluid pressure p, traction t, vector displacement u, and nonlinear/free (NL1 and NL2) boundary conditions considered in the analysis
In our previous study [11], we explained the treatment of the two free boundary problems that arise in the model, which we briefly describe here. First, a traction equal to the fluid pressure and a constant volumetric flow rate should be applied to the surface of the infusion cavity, however, the fluid pressure and cavity deformation are not known a priori. Second, at the beginning of the infusion the tissue is in contact with the catheter, but as the infusion proceeds the tissue near the infusion cavity begins to displace away from the catheter and axial flow develops in the gap between the catheter and the tissue. The model includes two specially formulated layers of elements at the tip and sides of the catheter to represent these conditions (Fig. 2).
Med Biol Eng Comput (2014) 52:841–849
843
The first layer (NL1) allows for the modeling of a constant infusion flow rate, whereas the other layer (NL2) allows for modeling the Poiseuille flow in the annular region [11]. The remaining boundary conditions are illustrated in Fig. 2. For the consideration of the constraint caused by the material surrounding the computational domain, a fixed radial displacement was assumed on the outer surface of the cylinder and a fixed condition on the external surface of the hemisphere, which are different from the zero-traction conditions on these boundaries assumed in our previous study [11]. To assess the influence of this difference, a comparative analysis was performed considering a tractionfree external boundary.
It should be noted that the nonlinear parameters αi describe the degree of nonlinearity of the stress–strain curve and, thus, the change of the shear modulus with strain. They also indirectly account for variations of the shear modulus due to the changes in the void ratio with strain. Consistent with the data reported by Miller and Chinzei [17], one term was considered in the series, with α1 taken to be −4.71. The Poisson’s ratio was taken to be 0.35, in agreement with previous analyses [7, 12, 25–27]. Consistent with analyses of infusion into brain tissue or brain tissue phantoms [7, 11, 12, 14, 25–27], the hydraulic conductivity was assumed to depend on tissue dilatation e as
2.2 Governing equations
κ = κ0 exp(M e),
Brain tissue was represented as a biphasic medium consisting of solid and fluid phases. This theoretical framework assumes that both phases are intrinsically incompressible but that the medium may compress by expulsion of the fluid. The governing equations are force equilibrium and mass conservation of the mixture [1, 11, 27], which may be, respectively expressed as
∇ · (se − pI) = 0
(2)
∇ · (vs − κ∇p) = 0,
(3)
where se is effective Cauchy stress tensor, p is the interstitial fluid pressure, I is the identity tensor, vs is the velocity of the solid matrix, and κ is the hydraulic conductivity. Whereas a linear elastic model is suitable for modeling backflow in agarose gel [11], experimental studies have demonstrated that the solid phase of brain tissue exhibits nonlinear stretch-stress curves and have described such behavior using an Ogden-type material model [10, 15–17]. Thus, we represented the solid phase of the infused domain with the Ogden-type hyperelastic function
W=
N 2µi αi 1 + α2 i + α3 i − 3 + (J −αi βi − 1)/βi , (4) 2 αi i=1
where λ1, λ2, λ3 are the principal stretch ratios, αi, μi and βi are material parameters, and J is the determinant of the deformation gradient tensor. The coefficients μi are related to the initial shear modulus G0 as
G0 =
N
µi ,
(5)
i=1
and the parameters βi are related to the Poisson’s ratio v by
βi =
ν . 1 − 2ν
(6)
(7)
where κ0 is the hydraulic conductivity at zero strain and M is a non-dimensional parameter. The elements that comprise the boundary layer NL2 that represents the annular region that develops outside the catheter were represented using the same energy function of Eq. (4) using one term of the series. The Poisson’s ratio was assumed to be zero to allow a high degree of compressibility. The other two parameters of the function, namely the nonlinear parameter α(BF) and the shear modulus G(BF), were determined by fitting experimental data from agarose gels, as explained in the next subsection. 2.3 Model calibration and validation The FE model was calibrated with in vitro data from agarose gel experiments reported by Raghavan et al. [20] since there is limited in vivo data from experimental or clinical infusion studies. First, the fixed values of the shear modulus G(BF) and the parameter α(BF) of the boundary layer NL2 were manually varied so that the length of the backflow zone and the fluid pressure from the model were fitted with the experimental data for a 0.98-mm radius catheter [20] by minimizing the objective function error f defined as
n � f = i=1
�
piEXP − piFEM piEXP
�2
+
n � i=1
�
i i − LFEM LEXP i LEXP
�2 1/2
,
(8)
i i where pEXP and pFEM are the experimental and numerical i i fluid pressure values, respectively, while LEXP and LFEM are the experimental and numerical backflow lengths, respectively. Since nonlinear material properties have not been reported for agarose gels, this calibration was performed representing the agarose gel as an elastic material using an initial shear modulus G0 of 2,300 Pa and a Poisson’s ratio equal to 0.45, consistent with our previous study [11]. Since the end of the backflow zone from the FE model is i not clearly defined, the backflow length LFEM was defined
13
844
Med Biol Eng Comput (2014) 52:841–849
to be the distance along the catheter from its tip to the point where the fluid velocity was 3 % of the maximum fluid velocity, as this percentage provided the best fit of the experimental data of Raghavan et al. [20]. This is different from our previous study [11], which used a displacement criterion, since the former criterion yielded backflow lengths much smaller than those predicted by the contours of the fluid pressure and fluid velocity. After calibrating the backflow-layer parameters with the experimental measurements for the 0.98-mm radius catheter, validation of the model was performed by comparing FE predictions with the experimental results of Raghavan et al. [20] for 0.33-mm radius and 0.50-mm radius catheters. In addition to predicting the backflow zone, this model can describe infusions when the fluid distribution is approximately spherical. For these cases, the reported backflow distance quantifies the length of the velocity distribution from the catheter tip toward the surface of the tissue along the lateral surface of the catheter. 2.4 Fitting of a power‑law formula The calibrated and validated FE model was used to fit the following simple power-law formula for the backflow length L
L = c Qw κ0x rcy Gz0 ,
(9)
which is similar to the formula [Eq. (1)] developed by Morrison et al. [18]. In this power-law formula, the parameters c, w, x, y, and z were obtained by fitting results of the FE model using the minimization function fmin of the program MATLAB (MathWorks Inc., Natick, MA). Note that the hydraulic resistance R in Eq. (1) is inversely proportional to the hydraulic conductivity κ. We chose to include the initial hydraulic conductivity κ0 directly in the scaling relationship in Eq. (9), as the hydraulic conductivity κ varies with deformation in our nonlinear model. In addition, since the infusion fluid is typically physiological saline, we do not expect there to be significant variation in the fluid viscosity µ and, therefore, it was omitted it from the scaling relationship. The data to fit the power-law formula was obtained by running the FE model for combinations of parameters within the ranges shown in Table 1. For the initial shear modulus G0, two ranges were chosen, where the first range
covers the values that have been reported for brain tissue [10, 15–17, 24, 29] and the second range considers a stiffer material in order to yield infinitesimal deformations more consistent with the closed-form solution. The latter range will allow for comparisons with other studies that have assumed that the brain tissue is completely rigid. For the sake of clarity in the remaining presentation, the formula of Eq. (1) developed by Morrison et al. [18] will be referred as the closed-form power-law formula and the formula of Eq. (9) fitted with results of the FE model will be described as the numerical power-law formula. 2.5 Evaluation of the numerical power‑law formula To determine the ability of numerical power-law formula to estimate the FE results, predictions of the numerical formula were compared with independent FE results for 21 random combinations of parameters chosen within the ranges of Table 1 for both the stiffer and compliant materials. First, a coefficient of determination (R2) was calculated to measure the association between the numerical powerlaw formula and these FE results. In addition, a Bland– Altman analysis [2] was performed to assess how well the numerical power-law formula predicts the FE results. 2.6 Assessment of effect of geometric and material nonlinearities Since the relaxed shear modulus of brain tissue has been documented to be less than 1,000 Pa [15–17], to better understand the influence of geometric nonlinearities associated with large deformations around the catheter tip, another set of 30 simulations were performed by varying G0 within the range 200 to 1,000 Pa, while the other parameters were taken near the average value of the range used to fit the power-law formulas (that is, rc = 0.5 mm, Q = 4.5 µL min−1, and κ0 = 4.0 mm4 N−1 s−1). To assess the influence of material nonlinearities, analyses were also performed by setting the parameter α1 of Eq. (4) equal to 1, as this parameter modifies the degree of nonlinearity of the function and a value equal to 1 describes a more linear relation between stretch and strain. Finally, to analyze the influence of the nonlinear variation of hydraulic conductivity with strain, these simulations were also performed by setting the parameter M of Eq. (7) equal to zero.
Table 1 Range of parameters used in the FE simulations that were used to fit the numerical power-law formula [Eq. (9)] Parameter
Number of samples
G0 (Pa)
Q (µL min−1)
κ0 (mm4 N−1 s−1)
rc (mm)
Brain G0
281
200–4,000
0.2–8
2–6
0.105, 0.33, 0.50, 0.98
High G0
201
10,000–200,000
0.2–8
2–6
0.105, 0.33, 0.50, 0.98
13
Med Biol Eng Comput (2014) 52:841–849
3 Results
845
3.2 Fitting and evaluation of the numerical power‑law formula
3.1 Model validation For the 0.98-mm radius catheter and flow rates in the range 1–8 µL min−1, the simulations with a zero-traction condition on the outer boundary [11] yielded maximum differences of 3.0 and 1.6 % for the backflow length and maximum fluid velocity, respectively, with respect to the simulations with the zero-displacement boundary conditions adopted in this study, which represent the restraint caused by the material adjacent to the computational domain. The best fitting of the experimental data (Fig. 3) for agarose gels [20] was obtained with a fixed shear modulus G(BL) of the boundary layer equal to 120 Pa and a fixed nonlinear parameter α(BL) equal to 2.3 (R2 = 0.981 and 0.995 for the fluid pressure and backflow length, respectively). With this set of constants, the FE model predicted well the independent experimental data for agarose gels [20] for both the 0.33-mm radius catheter (R2 = 0.996 and 0.999 for the fluid pressure and backflow length, respectively) and 0.50-mm radius catheter (R2 = 0.988 and 0.981 for the fluid pressure and backflow length, respectively). In contrast to the fitting obtained by Raghavan et al. [20, Table 3] in which a different shear modulus was needed for each catheter radius, this model allowed predicting those experimental results with a single set of material properties.
For the FE simulations with the higher range of the initial shear modulus G0, the numerical formula for backflow length L was
L = 1.50 Q0.76 κ0−0.99 rc0.43 G−0.60 , 0
(10)
whereas the numerical formula fitted with the FE simulations for the lower range of the initial shear modulus G0 was
L = 19.76 Q0.62 κ0−0.87 rc0.13 G−0.078 , 0
(11)
where L is given in mm. Larger differences between the numerical and closed-form solution exponents were obtained for the lower range of G0 (Table 2), as the exponents that affect the catheter radius and initial shear modulus were only 15 and 12 % of the closed-form values, respectively. The predictions of backflow length obtained with the numerical power-law formula were closer to the FE results for stiffer materials (10,000 Pa
ORIGINAL ARTICLE
Backflow length predictions during flow‑controlled infusions using a nonlinear biphasic finite element model Gustavo A. Orozco · Joshua H. Smith · José J. García
Received: 20 November 2013 / Accepted: 13 August 2014 / Published online: 26 August 2014 © International Federation for Medical and Biological Engineering 2014
Abstract A previously proposed finite element model that considers geometric and material nonlinearities and the free boundary problems that occur at the catheter tip and in the annular zone around the lateral surface of the catheter was revised and was used to fit a power-law formula to predict backflow length during infusions into brain tissue. Compared to a closed-form solution based on linear elasticity, the power-law formula for compliant materials predicted a substantial lower influence of the shear modulus and catheter radius on the backflow length, whereas the corresponding influence for stiffer materials was more consistent with the closed-form solution. The finite element model predicted decreases of the backflow length for reduction of the shear modulus for highly compliant materials (shear modulus less than 500 Pa) due to the increased area of infusion and the high fluid fraction near the infusion cavity that greatly increased the surface area available for fluid transfer and reduced the hydraulic resistance toward the tissue. These results show the importance of taking into account the material and geometrical nonlinearities that arise near the infusion surface as well as the change of hydraulic conductivity with strain for a proper G. A. Orozco Escuela de Ingeniería Mecánica, Universidad del Valle, Cali, Colombia e-mail: [email protected] J. H. Smith (*) Department of Mechanical Engineering, Lafayette College, Easton, PA 18042, USA e-mail: [email protected] J. J. García Escuela de Ingeniería Civil y Geomática, Universidad del Valle, Cali, Colombia e-mail: [email protected]
characterization of backflow length during flow-controlled infusions into the brain. Keywords Convection-enhanced delivery · CED · Infusion · Biphasic materials
1 Introduction Convection-enhanced delivery (CED) is a means to bypass the blood–brain barrier and deliver controlled volumes of drugs over localized zones of the brain for the treatment of diseases and tumors [3, 19]. Many experimental studies have been performed in agarose gels and animals to understand the effect of different parameters on the outcome of CED. Early brain infusion experiments with animals at higher flow rates showed backflow, in which an annular zone is formed outside the catheter and the infused drug preferentially flows toward the surface of the brain rather than through the tissue toward the area targeted for delivery [6]. In this case, the transport of the infused drug mainly occurs from the annular zone around the catheter rather than from its tip. To provide recommendations for clinical studies, controlled experiments using agarose gels as a brain tissue phantom have been performed for the measurement of infusion pressure, drug distribution, and backflow length [8, 9]. The foregoing experimental and clinical studies have been complemented with analytical models aimed to study the variables that affect backflow. The seminal model by Morrison et al. [18], later enhanced by Raghavan et al. [20], allows for the prediction of the backflow length and fluid pressure by considering the deformation of the tissue around the external boundary of the catheter, the axial flow in the annular gap developed around the cannula, and the
13
842
Med Biol Eng Comput (2014) 52:841–849
radial flow from this annular region into the porous tissue. Specifically, Morrison et al. [18] developed a power-law correlation for the backflow length
L ∝ Q0.6 R0.8 rc0.8 G−0.6 µ−0.2 , 0
(1)
where Q is the infusion flow rate, R is a tissue hydraulic resistance, rc is the catheter radius, G0 is the tissue shear modulus, and µ is the fluid viscosity. However, this formula was derived under some limiting assumptions, such as considering the solid phase of the infused tissue as a linearly elastic material under infinitesimal deformations, whereas mechanical testing has shown large deformations under physiological loadings [10, 15–17]. As recommendations for defining the parameters of infusions during clinical studies are being influenced by the backflow predictions of the analytical models [21–23] and as experimental infusion studies in animals have shown that deformation of the brain may be significant [28], it is important to develop models that incorporate the nonlinear behavior of brain tissue. To relax the assumptions of the analytical solution and to include geometric and material nonlinearities documented in experiments, we developed a finite element (FE) model of backflow in brain tissue phantoms [11] that was validated with published experimental results [20]. In this study, we updated our finite element model of backflow with a material model for brain tissue that is capable of reproducing the nonlinear stress–strain behavior that has been documented experimentally [10, 15–17]. We then used our revised finite element model to develop a powerlaw type formula similar to Eq. (1) that relates backflow length to the infusion flow rate and properties of the infused tissue. Lastly, we assessed the influence of nonlinear effects in the development of backflow.
2 Methods 2.1 Geometry and boundary conditions Our FE model of infusion into brain tissue from catheters of 0.105, 0.33, 0.50, and 0.98-mm radius was developed using ABAQUS 6.9-1 (Simulia, Providence, RI) considering axial symmetry and the planar geometry of Fig. 1. The computational domain extended between the outer radius of the catheter and a cylinder of 20-mm radius, with the lower portion of the model represented as a 20-mm radius hemisphere, consistent with previous studies [11, 18, 20]. Following a mesh convergence study based on the maximum fluid pressure, the domain was meshed with 4,486, 3,431, 3,364, and 3,241 elements of type CAX4P for the 0.105, 0.33, 0.50, and 0.98-mm radius catheters, respectively. These elements form an axisymmetric domain with respect to the vertical z-axis.
13
Fig. 1 Sketch of the domain considered in the FE model of backflow around a 0.98-mm radius catheter
Fig. 2 Representation of the radial displacement ur, radial fluid velocity qr, fluid pressure p, traction t, vector displacement u, and nonlinear/free (NL1 and NL2) boundary conditions considered in the analysis
In our previous study [11], we explained the treatment of the two free boundary problems that arise in the model, which we briefly describe here. First, a traction equal to the fluid pressure and a constant volumetric flow rate should be applied to the surface of the infusion cavity, however, the fluid pressure and cavity deformation are not known a priori. Second, at the beginning of the infusion the tissue is in contact with the catheter, but as the infusion proceeds the tissue near the infusion cavity begins to displace away from the catheter and axial flow develops in the gap between the catheter and the tissue. The model includes two specially formulated layers of elements at the tip and sides of the catheter to represent these conditions (Fig. 2).
Med Biol Eng Comput (2014) 52:841–849
843
The first layer (NL1) allows for the modeling of a constant infusion flow rate, whereas the other layer (NL2) allows for modeling the Poiseuille flow in the annular region [11]. The remaining boundary conditions are illustrated in Fig. 2. For the consideration of the constraint caused by the material surrounding the computational domain, a fixed radial displacement was assumed on the outer surface of the cylinder and a fixed condition on the external surface of the hemisphere, which are different from the zero-traction conditions on these boundaries assumed in our previous study [11]. To assess the influence of this difference, a comparative analysis was performed considering a tractionfree external boundary.
It should be noted that the nonlinear parameters αi describe the degree of nonlinearity of the stress–strain curve and, thus, the change of the shear modulus with strain. They also indirectly account for variations of the shear modulus due to the changes in the void ratio with strain. Consistent with the data reported by Miller and Chinzei [17], one term was considered in the series, with α1 taken to be −4.71. The Poisson’s ratio was taken to be 0.35, in agreement with previous analyses [7, 12, 25–27]. Consistent with analyses of infusion into brain tissue or brain tissue phantoms [7, 11, 12, 14, 25–27], the hydraulic conductivity was assumed to depend on tissue dilatation e as
2.2 Governing equations
κ = κ0 exp(M e),
Brain tissue was represented as a biphasic medium consisting of solid and fluid phases. This theoretical framework assumes that both phases are intrinsically incompressible but that the medium may compress by expulsion of the fluid. The governing equations are force equilibrium and mass conservation of the mixture [1, 11, 27], which may be, respectively expressed as
∇ · (se − pI) = 0
(2)
∇ · (vs − κ∇p) = 0,
(3)
where se is effective Cauchy stress tensor, p is the interstitial fluid pressure, I is the identity tensor, vs is the velocity of the solid matrix, and κ is the hydraulic conductivity. Whereas a linear elastic model is suitable for modeling backflow in agarose gel [11], experimental studies have demonstrated that the solid phase of brain tissue exhibits nonlinear stretch-stress curves and have described such behavior using an Ogden-type material model [10, 15–17]. Thus, we represented the solid phase of the infused domain with the Ogden-type hyperelastic function
W=
N 2µi αi 1 + α2 i + α3 i − 3 + (J −αi βi − 1)/βi , (4) 2 αi i=1
where λ1, λ2, λ3 are the principal stretch ratios, αi, μi and βi are material parameters, and J is the determinant of the deformation gradient tensor. The coefficients μi are related to the initial shear modulus G0 as
G0 =
N
µi ,
(5)
i=1
and the parameters βi are related to the Poisson’s ratio v by
βi =
ν . 1 − 2ν
(6)
(7)
where κ0 is the hydraulic conductivity at zero strain and M is a non-dimensional parameter. The elements that comprise the boundary layer NL2 that represents the annular region that develops outside the catheter were represented using the same energy function of Eq. (4) using one term of the series. The Poisson’s ratio was assumed to be zero to allow a high degree of compressibility. The other two parameters of the function, namely the nonlinear parameter α(BF) and the shear modulus G(BF), were determined by fitting experimental data from agarose gels, as explained in the next subsection. 2.3 Model calibration and validation The FE model was calibrated with in vitro data from agarose gel experiments reported by Raghavan et al. [20] since there is limited in vivo data from experimental or clinical infusion studies. First, the fixed values of the shear modulus G(BF) and the parameter α(BF) of the boundary layer NL2 were manually varied so that the length of the backflow zone and the fluid pressure from the model were fitted with the experimental data for a 0.98-mm radius catheter [20] by minimizing the objective function error f defined as
n � f = i=1
�
piEXP − piFEM piEXP
�2
+
n � i=1
�
i i − LFEM LEXP i LEXP
�2 1/2
,
(8)
i i where pEXP and pFEM are the experimental and numerical i i fluid pressure values, respectively, while LEXP and LFEM are the experimental and numerical backflow lengths, respectively. Since nonlinear material properties have not been reported for agarose gels, this calibration was performed representing the agarose gel as an elastic material using an initial shear modulus G0 of 2,300 Pa and a Poisson’s ratio equal to 0.45, consistent with our previous study [11]. Since the end of the backflow zone from the FE model is i not clearly defined, the backflow length LFEM was defined
13
844
Med Biol Eng Comput (2014) 52:841–849
to be the distance along the catheter from its tip to the point where the fluid velocity was 3 % of the maximum fluid velocity, as this percentage provided the best fit of the experimental data of Raghavan et al. [20]. This is different from our previous study [11], which used a displacement criterion, since the former criterion yielded backflow lengths much smaller than those predicted by the contours of the fluid pressure and fluid velocity. After calibrating the backflow-layer parameters with the experimental measurements for the 0.98-mm radius catheter, validation of the model was performed by comparing FE predictions with the experimental results of Raghavan et al. [20] for 0.33-mm radius and 0.50-mm radius catheters. In addition to predicting the backflow zone, this model can describe infusions when the fluid distribution is approximately spherical. For these cases, the reported backflow distance quantifies the length of the velocity distribution from the catheter tip toward the surface of the tissue along the lateral surface of the catheter. 2.4 Fitting of a power‑law formula The calibrated and validated FE model was used to fit the following simple power-law formula for the backflow length L
L = c Qw κ0x rcy Gz0 ,
(9)
which is similar to the formula [Eq. (1)] developed by Morrison et al. [18]. In this power-law formula, the parameters c, w, x, y, and z were obtained by fitting results of the FE model using the minimization function fmin of the program MATLAB (MathWorks Inc., Natick, MA). Note that the hydraulic resistance R in Eq. (1) is inversely proportional to the hydraulic conductivity κ. We chose to include the initial hydraulic conductivity κ0 directly in the scaling relationship in Eq. (9), as the hydraulic conductivity κ varies with deformation in our nonlinear model. In addition, since the infusion fluid is typically physiological saline, we do not expect there to be significant variation in the fluid viscosity µ and, therefore, it was omitted it from the scaling relationship. The data to fit the power-law formula was obtained by running the FE model for combinations of parameters within the ranges shown in Table 1. For the initial shear modulus G0, two ranges were chosen, where the first range
covers the values that have been reported for brain tissue [10, 15–17, 24, 29] and the second range considers a stiffer material in order to yield infinitesimal deformations more consistent with the closed-form solution. The latter range will allow for comparisons with other studies that have assumed that the brain tissue is completely rigid. For the sake of clarity in the remaining presentation, the formula of Eq. (1) developed by Morrison et al. [18] will be referred as the closed-form power-law formula and the formula of Eq. (9) fitted with results of the FE model will be described as the numerical power-law formula. 2.5 Evaluation of the numerical power‑law formula To determine the ability of numerical power-law formula to estimate the FE results, predictions of the numerical formula were compared with independent FE results for 21 random combinations of parameters chosen within the ranges of Table 1 for both the stiffer and compliant materials. First, a coefficient of determination (R2) was calculated to measure the association between the numerical powerlaw formula and these FE results. In addition, a Bland– Altman analysis [2] was performed to assess how well the numerical power-law formula predicts the FE results. 2.6 Assessment of effect of geometric and material nonlinearities Since the relaxed shear modulus of brain tissue has been documented to be less than 1,000 Pa [15–17], to better understand the influence of geometric nonlinearities associated with large deformations around the catheter tip, another set of 30 simulations were performed by varying G0 within the range 200 to 1,000 Pa, while the other parameters were taken near the average value of the range used to fit the power-law formulas (that is, rc = 0.5 mm, Q = 4.5 µL min−1, and κ0 = 4.0 mm4 N−1 s−1). To assess the influence of material nonlinearities, analyses were also performed by setting the parameter α1 of Eq. (4) equal to 1, as this parameter modifies the degree of nonlinearity of the function and a value equal to 1 describes a more linear relation between stretch and strain. Finally, to analyze the influence of the nonlinear variation of hydraulic conductivity with strain, these simulations were also performed by setting the parameter M of Eq. (7) equal to zero.
Table 1 Range of parameters used in the FE simulations that were used to fit the numerical power-law formula [Eq. (9)] Parameter
Number of samples
G0 (Pa)
Q (µL min−1)
κ0 (mm4 N−1 s−1)
rc (mm)
Brain G0
281
200–4,000
0.2–8
2–6
0.105, 0.33, 0.50, 0.98
High G0
201
10,000–200,000
0.2–8
2–6
0.105, 0.33, 0.50, 0.98
13
Med Biol Eng Comput (2014) 52:841–849
3 Results
845
3.2 Fitting and evaluation of the numerical power‑law formula
3.1 Model validation For the 0.98-mm radius catheter and flow rates in the range 1–8 µL min−1, the simulations with a zero-traction condition on the outer boundary [11] yielded maximum differences of 3.0 and 1.6 % for the backflow length and maximum fluid velocity, respectively, with respect to the simulations with the zero-displacement boundary conditions adopted in this study, which represent the restraint caused by the material adjacent to the computational domain. The best fitting of the experimental data (Fig. 3) for agarose gels [20] was obtained with a fixed shear modulus G(BL) of the boundary layer equal to 120 Pa and a fixed nonlinear parameter α(BL) equal to 2.3 (R2 = 0.981 and 0.995 for the fluid pressure and backflow length, respectively). With this set of constants, the FE model predicted well the independent experimental data for agarose gels [20] for both the 0.33-mm radius catheter (R2 = 0.996 and 0.999 for the fluid pressure and backflow length, respectively) and 0.50-mm radius catheter (R2 = 0.988 and 0.981 for the fluid pressure and backflow length, respectively). In contrast to the fitting obtained by Raghavan et al. [20, Table 3] in which a different shear modulus was needed for each catheter radius, this model allowed predicting those experimental results with a single set of material properties.
For the FE simulations with the higher range of the initial shear modulus G0, the numerical formula for backflow length L was
L = 1.50 Q0.76 κ0−0.99 rc0.43 G−0.60 , 0
(10)
whereas the numerical formula fitted with the FE simulations for the lower range of the initial shear modulus G0 was
L = 19.76 Q0.62 κ0−0.87 rc0.13 G−0.078 , 0
(11)
where L is given in mm. Larger differences between the numerical and closed-form solution exponents were obtained for the lower range of G0 (Table 2), as the exponents that affect the catheter radius and initial shear modulus were only 15 and 12 % of the closed-form values, respectively. The predictions of backflow length obtained with the numerical power-law formula were closer to the FE results for stiffer materials (10,000 Pa
