Optimization of intraperitoneal aerosolized drug delivery using computational fluid dynamics (CFD) modeling

Intraperitoneal (IP) aerosolized anticancer drug delivery was recently introduced in the treatment of patients with peritoneal metastases. However, little is known on the effect of treatment parameters on the spatial distribution of the aerosol droplets in the peritoneal cavity. Here, computational fluid dynamics (CFD) modeling was used in conjunction with experimental validation in order to investigate the effect of droplet size, liquid flow rate and viscosity, and the addition of an electrostatic field on the homogeneity of IP aerosol. We found that spatial distribution is optimal with small droplet sizes (1–5 µm). Using the current clinically used technology (droplet size of 30 µm), the optimal spatial distribution of aerosol is obtained with a liquid flow rate of 0.6 mL s−1. Compared to saline, nebulization of higher viscosity liquids results in less homogeneous aerosol distribution. The addition of electrostatic precipitation significantly improves homogeneity of aerosol distribution, but no further improvement is obtained with voltages higher than 6.5 kV. The results of the current study will allow to choose treatment parameters and settings in order to optimize spatial distribution of IP aerosolized drug, with a potential to enhance its anticancer effect.


Materials and methods
In-vitro experiments. In-vitro experiments were performed using a plexiglass box (185 × 135 × 152 mm 3 , Fig. 2), as previously described 17 . The nebulizer was placed at the top surface of the box through a GelPOINT Mini access platform (Applied Medical, Amersfoort, The Netherlands). A volume of 20 mL of black ink (Pelikan nv, Groot-Bijgaarden, Belgium) was nebulized with a volumetric flow rate (Q) of 0.5 mL s −1 and upstream pressure of 20 bar using a high pressure injector (Injektron 82 M, Medtron AG, Saarbrücken, Germany). Fresh rat peritoneal tissue was obtained from animals that were culled after neurosurgical experimentation at Ghent University. Tissue samples (approximately 2 mm thick and 2 by 2 cm) were positioned on metal plates located at four different sites in the box to allow visualizing and quantifying the distribution of black ink aerosol droplets:  . The samples were positioned with the peritoneal surface facing the aerosol, and the skin surface in contact with the metal plates. A pressure regulated insufflator was used to create a CO 2 pressure of 12 mmHg (1600 Pa) in the box prior to black ink nebulization. After each experiment (with a duration of 30 min), digital images were obtained to document ink coverage of the tissue samples. Subsequently, the extent of ink coverage and distribution were calculated by means of threshold functions (brightness 1-100) with ImageJ software (version 1.51, National Institutes of Health, Bethesda, Maryland, United States, available from https:// imagej. nih. gov/ ij/ index. html). A region of interest (ROI) was drawn around the tissue border and the area of the ROI was measured. The proportion (%) of stained tissue surfaces was calculated as follows: The measurements were performed in triplicate.
Characteristics of aerosol droplets. Aerosol droplet size measurement. Measurements of the volume weighed particle size distribution (PSD) of aerosol droplets was performed using laser diffraction (Mastersizer S long bench, Malvern Instruments, Malvern, United Kingdom). An open laser beam (water vs. air, refractive index of 1.33 and 1.00, respectively) was created with a 300F lens (0.5-900 μm) after 10 s of nebulization. The laser diffraction instrument obtains the size distribution of aerosol droplets by measuring the angular variation in intensity of light scattered as a laser beam passes through the droplets. When the aerosol droplets interfere with the laser beam, they create a diffraction pattern. Droplet size measurements were carried out at 3.5 and 10 cm from the tip of the nozzle and the lens, respectively. Several experiments were performed to cover a range of flow rates (0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 mL s −1 ). Results were exported as the median of the volume distribution, D(v,0.5), i.e. the volume median diameter at which 50 vol% of the aerosol droplets were either finer or coarser than the predicted value, with standard deviation.
Cone angle of nebulization. The nozzle geometry and characteristics affect the aerosol droplet behavior during nebulization 18 . The diameter of the orifice and the cone angle of nebulization are two important characteristics of the nozzle. The cone angle has a relationship with the driving pressure, flow rate, and viscosity of the fluid. In the current PIPAC technique, the high-pressure nozzle has an orifice with a diameter of 200 µm. The experiments were performed using water and flow rates of 0.2, 0.5, 0.8 and 1.1 mL s −1 with a maximal pressure of 20 bar. During each experiment, digital images were obtained during nebulization and imported to ImageJ to measure the cone angle. During CFD simulation, the cone angle was determined directly by COMSOL Multiphysics (Measure accumulator) by defining two edges drawn in a plane from the nozzle tip to the outer periphery of the aerosol.
Viscosity measurement. In order to investigate the effect of the viscosity of the carrier liquid, we measured the viscosity of Icodextrin (Baxter Healthcare Ltd, Illinois, US), a glucose polymer preparation, at concentrations of 4% and 7.5% using a capillary viscometer (Paragon Scientific Ltd, Birkenhead, UK). The procedure was performed according to the guidance from the European Pharmacopoeia 10.0 19 . Saline solution (NaCl 0.9%) was used to calibrate of the viscometer. The experiments were performed in triplicate.
Electrostatic precipitation. To model electrostatic precipitation, we used a high voltage power supply (LD Didactic, Hürth, Germany, 0-10 kV, current 2 mA). An active cable with a steel brush electrode and four return electrodes connected to four metal grounding plates were added to the box model (Fig. 2). The electrostatic precipitation was initiated from the start of nebulization.
Computational approach. Geometry and mesh generation of the box model. A 3D model of a rectangular box (185 × 135 × 152 mm 3 ) (Fig. 3A) was created using COMSOL Multiphysics (COMSOL, Inc., Burlington, VT, USA). The fluid domain within the box geometry was spatially discretized using tri/tetrahedral grids with local refinements near the nozzle inlet and locations A-D (Fig. 3B). To ensure valid results, a mesh sensitivity study was performed, resulting in an optimal volume mesh containing 98,798 tetrahedral elements. The physical domain was discretized into a large number of tri/tetrahedral computational elements with local refinements near critical locations, i.e. near the nozzle inlet and near the tissue holding plates inside the box. The mesh quality was evaluated by examining the skewness measure of mesh elements which is based on equiangular skew. It was assumed that a skewness is less than 0.9 was acceptable. The maximum skewness ranged from 0.68 to 0.85, in range with published recommendations 20 . Elements with a quality below 0.1 were considered as a poor quality mesh. To ensure valid results, a mesh sensitivity study was performed, resulting in an optimal hybrid volume mesh containing 98,798 tri/tetrahedral elements. The model was tested for three different grid densities, i.e. containing 48,088 (Mesh 1), 98,798 (Mesh 2) and 162,820 (Mesh 3) volume elements, by comparing the CO 2 phase stability after reaching a pressure of 12 mm Hg. Three cut-lines were considered in the X, Y and Z directions inside the box model. The average velocity of the CO 2 phase in these three directions inside the cavity was selected as one of the criteria for mesh independency. The comparison of the average velocity of CO2 phase for all three grid resolutions is provided in Fig. 4A. Comparing the results showed a 5.6% difference between Mesh 1 and Mesh 2, and a 0.7% difference between Mesh 2 and Mesh 3. The increase of element number from 98,798 to 162,820 grids did not change the results. Additionally, to investigate the effects of mesh resolution on aerosol   www.nature.com/scientificreports/ droplet distribution, the droplet deposition efficiency of 2 mesh densities, 98,798 (Mesh 2) and 162,820 (Mesh 3) elements, were compared for 4 regions of the in vitro box model (regions 1-4) for two different droplet diameters (1 and 30 µm). Comparison of the droplet deposition percentage over specific regions (region 1, 2, 3, and 4) of two mesh densities showed that the maximum error was 0.3% (Fig. 4B). Therefore, we adopted the mesh that contained 98,798 grid cells, the combination that gave the best grid independency and stability in trade-off with the computational cost. The CFD simulations were performed on the high performance computing infrastructure of Ghent University, and on a Dell Precision 5820 workstation with 256 GB RAM and 4.00 GHZ CPU, resulting in typical calculation times of 48-120 h.
Continuous phase modeling (CO 2 gas). As a first step, the CFD approach filled the virtual box with CO 2 gas to a pressure of 12 mmHg (1600 Pa). The governing equations of the CO 2 phase were considered under unsteady, laminar, incompressible flow and Newtonian fluid assumptions. The associated continuity and momentum (Navier-Stokes) equations can be expressed as follows using Einstein notation 21 : • Continuity equation  22 . Impacted droplets were assumed to adhere to the sidewalls if the distance between the droplet center and the sidewall was less than the droplet's diameter. Two-phase model. The simulation entailed a discrete (aerosol droplets) and continuous (gas) phase. Aerosol transport was modelled by tracking individual droplets as they move through the CO 2 gas using a Lagrangian approach 23 . Based on the volume fraction of the discrete phase, the two-phase flow can be described using either one-way or two-way coupling. The volume fraction is given by 24,25 : where N, v and V are the number of droplets, droplet volume, and the volume of gas (continuous phase) respectively. In the present study, the volume fraction ( ϕ ) is less than 10 -6 , indicating the continuous phase (CO 2 gas) significantly outweighs the discrete phase (aerosol droplets). Using this approach, the two-phase flow can be considered as dilute and the momentum imparted onto the fluid by the droplets can be neglected, justifying the assumption to use one-way coupled simulations 26,27 .
The simulation was solved in two stages. First, the CO 2 gas flow (continuous phase) was solved using a stationary study step. Then, the aerosol droplet (discrete phase) trajectories were computed using a time dependent study step. The solution from the stationary study was used to define the CO 2 pressure and velocity for the purpose of exerting a drag force on the droplets. The presence of droplets was not considered when modelling the CO 2 filling of the cavity.
Forces acting on individual aerosol droplets. Newton's second law describes the aerosol droplet transport in the Lagrangian formulation using a Cartesian coordinate system (for a given direction) as follows 28 : When adding electrostatic precipitation, an additional electrical force (F E ) is generated by the electrostatic field, which is proportional to the electrical charge of the aerosol droplets (q [C]) and the strength of the electrical field (E [V m −1 ]): The overall force balance (Eq. 4) can be rewritten in Cartesian coordinates as follows: Stokes and Weber numbers. When a volume of fluid is forced through a nozzle orifice, the pressure head of the fluid is converted into kinetic energy. The probability of deposition by inertial impaction is higher when the droplets are more likely to travel longer distances, which is based on the velocity and the size of the droplet. The Stokes number (Stk) is the ratio of the aerosol droplet stopping distance ( τ p = u − u p ) and half of the characteristic length (0.5 d N ) as follows 11,29 : where t f , d N  The Weber number (We) is the ratio between the inertial and the surface tension forces of liquid 30 : where σ [N m −1 ] is the surface tension. The We number indicates whether the kinetic or the surface tension energy is dominant: the higher the We number, the more dominant is the kinetic energy 30 . The value of the We number determines the mode of droplet breakup, reflecting the minimal initial inertial force required to cause droplet breakup assuming a certain restoring force, under an impulsive acceleration. The We number characterizes the tendency of the liquid to break up under competing inertia and surface tension forces.
In silico simulation. For CFD simulation of aerosol droplet transport, the same steps were applied in silico: insufflation of CO 2 to reach a pressure of 12 mmHg, followed by nebulization of ink. The computational domain was assumed to be incompressible. For the in vitro validation experiment and simulation of the box model, the internal pressure was increased by injecting CO 2 to resemble the time point at which the cavity is fully inflated at a steady state pressure of 12 mmHg. In this situation, the peritoneal cavity can be assumed to be an incompressible domain. Then, a drug-containing solution was aerosolized within the (inflated) peritoneal cavity by means of an atomizer. The CFD and Particle Tracing Modules were applied under appropriate initial and boundary conditions. The dimensions of the simulated tissue samples at locations A-D (see Fig. 3A) were 20 × 20 × 2 mm 3 . The locations of tissues in the in-vitro experiment and in-silico identical (A, B, C and D). Tissue properties (conductivity, density, Young's modulus, relative permeability and electrical conductivity) were defined using the material library model for human tissue in COMSOL Multiphysics. In order to estimate the deposition of aerosol on the simulated tissue samples, the following procedure was used to mimic the experimental (in vitro) setup: after running the simulation, the areas representing the virtual samples were selected and exported as a TIFF file. Using these files, . The simulations allowed to determine the effect of these parameters on the deposition of the aerosol droplets in different regions of the model. Condensation and evaporation effects were neglected for the PIPAC simulations. In pulmonary medicine, the elevated temperature and relative humidity of the airways are known to potentially affect the properties of aerosol droplets by causing evaporation, condensation, and hygroscopic growth 31 . However, using current laparoscopic technology, cold and dry CO 2 gas, which has a very low relative humidity (typically 0.0002%), is insufflated at room temperature (20 °C). Therefore, these two phenomenona were assumed to be negligible during PIPAC. The simulation was solved in two stages. First, the CO 2 gas flow (continuous phase) was solved for using a stationary study step. Then, the aerosol droplet (discrete phase) trajectories were computed using a time dependent study step. The solution from the stationary study was used to define the CO 2 pressure and velocity for the purpose of exerting a drag force on the droplets. The presence of droplets was not considered when modelling the CO 2 filling the cavity. A second-order (quadratic) discretization of fluid flow was used and the partial derivative equations were solved with a dominant second derivative term. A fully coupled method was applied to generate the single set of algebraic equations for all the involved physical models, and implemented in a single iteration scheme which was repeated until convergence was reached. The model was considered converged when the estimated error in the iterative solver was below 10 -6 .

Results and discussion
Experimental validation of the CFD model. In a first step, we modelled filling of the box with CO 2 gas to a pressure of 12 mmHg (1600 Pa), Fig. 5. Next, we simulated the spatial distribution of aerosol droplets after nebulization for PIPAC and ePIPAC (Fig. 6A,B). As expected, the simulation showed that most aerosol droplets were deposited at the bottom region of the box due to gravity and inertial impaction during PIPAC. Figure 6C,D display the droplet distribution after experiments of PIPAC and ePIPAC, respectively. As can be seen, the spatial distribution of droplets after ePIPAC is more homogenous than PIPAC in the box model Fig. 6E compares the percentage of tissue surface ink staining (regions A-D) obtained in vitro with the CFD simulation, showing an overall strong agreement between both. Tissue at the bottom surface was almost completely covered by black ink after PIPAC. Tissue on plate C was partially stained, while staining of tissue on plate D was very limited. By imposing an electrostatic field in the box to mimick ePIPAC, the aerosol distribution improved and was more homogenous. Both the experiment and simulation results show a significantly better aerosol droplet deposition at the top of the box model (plate D) when an electrical force is created.
Effect of aerosol droplet size. The aerosol droplet size is a main parameter affecting the type of forces that impact aerosol droplet transport 32 . As mentioned in "Stokes and Weber numbers" (Eq. 10), the Stk number, is a measure of the influence of the inertial effects during droplet transport. Thus, Stk as well as Re p have to be considered carefully to capture more accurately the underlying physics of droplet deposition. Understanding the relationship between the aerosol droplet size and forces is helpful to predict the behavior of the aerosol droplet. As shown in Eq. (5), the drag force has a direct relationship with the aerosol droplet diameter. Also, the gravitational force has a direct relationship with the aerosol droplet diameter to the 3rd power (see Eq. 7). When the size of the aerosol droplet raises, its mass and subsequently the gravitational force also significantly increase. Thus, for large aerosol droplets, drag and especially gravitational forces play a significant role 32 .  www.nature.com/scientificreports/ For the analysis of the effect of aerosol droplet size, the CFD modelled box was divided in four regions, corresponding to dorsal (region 1) to ventral (region 4) regions of the peritoneal cavity. Five different aerosol droplet diameters were modelled (see Fig. 7). The simulations showed that aerosol droplet deposition is significantly affected by their diameter. As shown in Fig. 7, the effect of gravity is quite significant for droplets in the range of 30-50 µm. For smaller droplets, gravity has only a small effect on the deposition pattern. Figure 8 shows the simulated deposition of aerosol droplets, according to diameter, in the four defined regions. With increasing diameter, most droplets deposited in the dorsal region (region 1). With a diameter of 30 µm, 60% of the aerosol droplets deposited in region 1, 30% in region 2, 10% in region 3, and none in region 4. With a diameter of 50 µm, 78% of aerosol droplets deposited in region 1. The results showed that a droplet diameter ranging between 1 and 5 µm resulted in the most homogenous spatial distribution, although coverage of the most ventral region (region 4) remained virtually absent with this flow rate, due to less perturbation of aerosols. More perturbation in the aerosol droplet flow results a better distribution. Figure 9A,B compare the deposition of the aerosol droplets with different flow rates in four regions of the box using CFD simulations. In these comparisons, the simulation was done for two different diameters, while all other properties except flow rate of the nebulization  www.nature.com/scientificreports/ were kept constant. For a constant droplet size, increases in flow rate led to a higher initial velocity, concomitant increases in inertial impaction and kinetic energy, and extensive deposition in region 1 (opposite the nebulizer). However, very low flow rates (< 0.4 mL s −1 ) also resulted in deposition of a large majority of the droplets in region 1. This might probably explained by the inability of the low flow rate and energy to break up the liquid into droplets, and by the dominant effect of gravity, especially for heavy and large droplets, leading to sedimentation 33 . Also, low flow rates tend to decrease the nebulization angle (see "Effect of liquid flow rate on nebulization cone angle"). Inertial impaction and gravitational force seem to have the most effect on the droplet behavior. For a smaller diameter (1 µm), homogeneity of aerosol deposition was found to be optimal with a flow rate of 0.7 mL s −1 . For diameter of 30 µm, the inertial impaction and gravitational force are higher, and a flow rate of 0.5 mL s −1 shows a better distribution of droplets in the box model. Obviously, for larger droplets, the large majority of droplets are deposited on the opposite site of the nebulizer. According to Fig. 9A, due to more complex aerosol flows, some aerosol droplets deposited at region 4 for flow rates of 0.6 mL s −1 and 0.7 mL s −1 . However, for larger droplets, this region remained unexposed. Figure 10 demonstrates the experimentally measured volume-weighted distributed density of the aerosol droplet diameters for different flow rates of the nebulization. Interestingly, analysis of different flow rates at a fixed maximal upstream injection pressure of 20 bar, showed that higher flow rates of liquid resulted in a decrease of the volume median diameter of the droplets. When increasing the flow rate, D(v,0.5) values were measured as 48 ± 2 µm (p < 0.0001), 40 ± 1 µm (p < 0.05), 35 ± 2 µm (p < 0.0001), 29 ± 2 µm (p < 0.0001), 28 ± 2 µm (p < 0.05) and 28 ± 3 µm (p < 0.05) for flow rates of 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 mL s −1 , respectively.  Table 1 shows impaction  *Q). Each bar represent the deposition percentage of the aerosol droplets in the box model. It shows preferential accumulation of the aerosol droplets in the dorsal region (region 1) due to inertial impaction and gravity, while the ventral region (region 4) remains unexposed. To obtain a more homogenous distribution of aerosol droplets, the flow rate and droplet size should be considered together. Both affect droplet behavior and relevant forces. The impaction parameter is just one of the factors that may affect the droplet distribution pattern. When changing the impaction parameter, the deposition pattern of the droplets will change, but for a more realistic prediction, all effective parameters (e.g. droplet size, flow rate, viscosity, surface tension and temperature gradient) should be considered. Here, we change only the values of the impaction parameters. For larger droplets, the extent of deposition of droplets in the dorsal region increases, due to inertial impaction and gravitational force. In clinical practice, it was recently recommended to increase the flow rate from 0.5 mL s −1 to 0.6 mL s −1 . On one hand, a higher flow rate of liquid will increase the initial velocity (assuming a constant droplet diameter), the inertial impaction and the deposition to the region opposite the nebulizer. On the other hand, increasing the flow rate from the nebulizer will lead to a smaller droplet size, which may improve spatial distribution of the aerosol. According to Fig. 10, the clinical device generates a mean droplet diameter range between 25 and 50 µm for different liquid flow rates. Results show that 0.6 mL s −1 can be an optimum flow rate of nebulization for the current PIPAC setup, considering the size of generated droplets. The homogeneity of distribution of aerosol droplets for 0.6 mL s −1 is slightly better than 0.5 mL s −1 . Consequently, the recommendation for increasing the flow rate to 0.6 mL s −1 may help to obtain a more homogenous distribution of aerosol droplets during PIPAC, although the effect is rather modest.

Effect of liquid flow rate. Fixed droplet diameter.
Effect of liquid flow rate on nebulization cone angle. Figure    www.nature.com/scientificreports/ measure accumulator was defined in the COMSOL Multiphysics and two edges were drawn from the nozzle tip to the outer periphery of the spray. As shown, at first, the angle increases significantly with increasing flow rates until reaching its maximum value after which it remains almost constant. Indeed, the cone angle for current PIPAC strategies with a maximum pressure of 20 bar and a flow rate of 0.5 mL s −1 (or 0.6 mL s −1 ) is around 70°. As can be seen in this figure, there is an exponential relationship between cone angle and flow rate as well.  Figures 13 and 14 show the results of CFD simulations of aerosol droplets being nebulized with different viscosities (saline (µ = 1 mPa s), Icodextrin 4% and 7.5%). As shown in these figures, at the highest viscosity (Icodextrin 7.5%), more aerosol droplets were deposited at the bottom (dorsal) region of the box model. The results of the CFD simulation for saline and Icodextrin solutions (4% and 7.5%) proved that the aerosol droplet distribution for liquid with lower viscosity is more homogenous than the liquid with higher viscosity.
Effect of electrical potential (voltage) during Epipac. Figure 15 displays the electrostatic field in the box model at different electrical potentials (4, 5, 6, 6.5, 7, 8 and 9 kV). The electrical field (E = v/m) is the ratio of electrical potential (voltage) and distance (meter) and consequently, electrical force has a direct relationship to   www.nature.com/scientificreports/ electrical potential. Increasing the electrical potential lead to an increase of the intensity of the electrostatic field in the box model. The charge of the droplets q were assumed to be − 1 for the simulation. It is obvious that the electrostatic force close to the brush electrode is higher than that further away. Figure 16 compares the in vitro experiments and CFD simulation results of the spatial distribution of black ink (the percentage of tissue surface ink staining for regions A-D) for different electrical potentials. Imposing an electrostatic field to the PIPAC setup (ePIPAC) results in a more homogenous aerosol droplet distribution in the box. A significant increase for proportion of black ink (a minimum of 60%) was observed at the stained tissue surface on the top wall of the box. The results showed an overall good agreement for the experiments and simulation. The electrical potential  www.nature.com/scientificreports/ sensitivity analysis revealed that the aerosol droplet distribution became more homogenous as the electrical potential increased, but no further improvements were obtained after 6.5 kV. Consequently, the optimum electrical potential for ePIPAC seemed to be 6.5 kV, considering black ink as nebulized liquid. Both the CFD model and the experimental results show a significantly better aerosol deposition at the top of the box (plate D) when the electrical potential reached 6.5 kV. It is notable that this effect depended on the type of aerosol droplets and their electrical charge.

Conclusions
Aerosolized intraperitoneal drug delivery holds considerable promise for the treatment of peritoneal cancer.
The results of the current study contribute to our understanding of the effects of different parameters on aerosol droplet behavior in the peritoneal cavity, and may guide a rational choice of treatment parameters during clinical PIPAC procedures. We found that spatial distribution is optimal with small droplet sizes (1-5 µm). Using the current clinically used technology (droplet size of 30 µm), the optimal spatial distribution of aerosol is obtained with a flow rate of 0.6 mL s −1 . The nebulization cone angle increases exponentially with flow rate, but a plateau is reached at 0.6 mL/s. Compared to saline, nebulization of higher viscosity liquids results in less homogeneous aerosol distribution. The addition of electrostatic precipitation significantly improves homogeneous aerosol distribution, but no further improvement is obtained with voltages higher than 6.5 kV. Further work will include the use of a realistic in vitro and CFD geometry, based on the actual anatomy, modelling of turbulence, and the use of high speed microscopic imaging to study the interaction of the aerosol droplets with the target tissues.