Optimization of Intraperitoneal Aerosolized Drug Delivery: A Computational Fluid Dynamics (CFD) and Experimental Study

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.


Introduction
Peritoneal metastases (PM) are a common manifestation of gastro-intestinal and gynecological cancers. Compared to other metastatic locations such as the liver, systemic chemotherapy is less active against PM, with a survival typically less than 10 months. [1] In addition, the quality of life of these patients is often poor due to debilitating symptoms such as obstruction or ascites formation. [2][3][4] Over the past decades, the role of intraperitoneal drug delivery (IPDD) has been explored as an addition to systemic chemotherapy. [5] During IPDD, PM are directly exposed to chemotherapy, while systemic toxic effects remain limited due to the pharmacological advantage conferred by the peritoneumplasma barrier. Different methods of IPDD are currently used, including catheter-based repeated intraperitoneal (IP) instillation and hyperthermic intraperitoneal chemoperfusion (HIPEC) in association with cytoreductive surgery (CRS). [6,7] However, many patients present with widespread and/or unresectable disease. Pressurized IntraPeritoneal Aerosol Chemotherapy (PIPAC) was recently introduced as an innovative drug delivery method to treat PM. [8,9] PIPAC combines IPDD with a minimally invasive approach (laparoscopy) and it is repeatable. Furthermore, PIPAC may improve drug penetration in tumoral tissue due to an increased intraperitoneal pressure. [10] During PIPAC, a CO2 pneumoperitoneum is first established (12 mmHg) to inflate the peritoneal cavity. Then, a drug-containing solution is typically aerosolized within the inflated peritoneal cavity by means of an atomizer (Capnopen ® , Capnomed, Zimmern, Germany) and a high-pressure injector (Injektron ™ 82 M, Medtron, Saarbrücken, Germany) with the aim to obtain a homogenous aerosol distribution within the peritoneal cavity. [11] (Figure 1)

Figure 1. Schematic illustration of pressurized intraperitoneal aerosol chemotherapy (PIPAC).
Several preclinical (in vitro [9] , ex vivo [12] , and in vivo [8] ) and clinical [13,14] studies demonstrated the feasibility, safety, and efficacy of PIPAC. However, the ideal procedural parameters such as flow rate, droplet size, and treatment duration are at present unknown. The aerosol droplet distribution is not uniform in the peritoneal cavity, probably due to the effects of inertial impaction and gravity. [15] The addition of an electrostatic force (electrostatic PIPAC or ePIPAC) may improve the spatial distribution of charged aerosol particles. [16] Here, we used computational fluid dynamics (CFD) modeling in combination with in vitro validation to define the optimal treatment parameters and thus guide clinical PIPAC therapy.

In-vitro experiments
In-vitro experiments were performed using a plexiglass box (185 x 135 x 152 mm 3 , Figure 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 samples 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: on the bottom surface without any additional obstacle (A), under a curved plastic obstacle on the bottom surface (B), on the side surface (left) (C), and on the top surface (D). A pressure regulated insufflator was used to create a CO2 pressure of 12 mmHg (1600 Pa) in the box prior to black ink nebulization. After each experiment (with a duration of 30 minutes), 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). 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: (%) = × 100

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 (Figure 2). The electrostatic precipitation was initiated from the start of nebulization.

Geometry and mesh generation of the box model
A 3D model of a rectangular box (185 x 135 x 152 mm³) ( Figure 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 ( Figure 3B). To ensure valid results, a mesh sensitivity study was performed, resulting in an optimal volume mesh containing 98,798 tetrahedral elements.
(A) (B) Figure 3. Illustration of the preparation of the simulation geometry for the CFD simulation approach: (A) 3D box model geometry, and (B) mesh generation.

Continuous phase modeling (CO2 gas)
As a first step, the CFD approach filled the virtual box with CO2 gas to a pressure of 12 mmHg (1600 Pa). The governing equations of the CO2 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 [20] : where ui [m s -1 ] (i = 1, 2, 3) are the flow velocity components, and xi [m] denotes the droplet position components.
 Momentum equation: is the viscous stress related to incompressible flow, ν [m 2 s -1 ] the fluid kinematic viscosity, p [Pa] the static pressure, and ρ0 [kg m -3 ] the fluid density. The inflow boundary condition for CO2 gas was assumed to be a pressure inlet (12 mmHg = 1600 Pa) at the level of the nozzle indicated on Figure 3. A no slip condition was used at the rigid walls. The Laminar Fluid Flow physics submodule within the CFD module in COMSOL Multiphysics™ was used to simulate the CO2 inflow (ρ= 1.977 kg m -3 ; ν= 7.44*10 -6 m 2 s -1 )

Discrete phase modeling (aerosol droplet)
After a stable pneumoperitoneum of 12 mmHg was created by means of the CO2 phase modelling, the second step of the CFD approach focused on analyzing aerosol droplet transport using the Particle Tracing for Fluid Flow module in COMSOL Multiphysics. Nebulization of a volume of 20 mL of black ink (density: 1070 kg m -³, dynamic viscosity: 2.1 mPa s) was simulated with a volumetric flow rate of 0.5 mL s -1 at a fixed injector position defined as an inlet at the middle of the top surface of the box. Aerosol droplets were injected during 40 seconds, and the total simulation time was 30 minutes. Based on the flow rate (Q=0.5 mL s -1 ) and the area of the orifice (AO) of the nebulizer nozzle (diameter = 200 µm), the initial velocity (u0) was calculated as 15.92 m s -1 (u0=Q/AO). The droplet diameter was set at 30 µm. [21] 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. Since the effect of the discrete phase on the continuous phase is negligible in the current model, a one-way coupling approach was used. [22] Aerosol transport was modelled by tracking individual droplets as they move through the CO2 gas using a Lagrangian approach [23] , which is based on the force balance on a single aerosol droplet.

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 [24] :  [N] are the drag, gravitational and electrical forces, respectively. FD is a resistance force related to the aerosol droplet characteristics (size and density) and the fluid viscosity. The resulting FD on a droplet with a diameter dp is defined by Stokes's equation as follows: since the Rep number affects the inertial impaction and drag force. The gravitational force (FG) due to gravitational acceleration (g [m² s -1 ]) can be expressed as follows: When adding electrostatic precipitation, an additional electrical force (FE) 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 (Equation 3) 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 ( = | − |) and half of the characteristic length (0.5 dN) as follows [11,25] : The Weber number (We) is the ratio between the inertial and the surface tension forces of liquid [26] : 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 [26] . 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.

Simulation
For CFD simulation of aerosol droplet transport, the same steps were applied in silico: insufflation of CO2 to reach a pressure of 12 mmHg, followed by nebulization of ink. 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 Figure 3A) were 20 x 20 x 2 mm³. A Measurement Accumulator was designed in the Particle Tracing Module to export the area of stained surface of the tissue. A threshold function was applied to determine the proportion of the surface of the tissue samples (A-D) covered by ink using ImageJ. The number of deposited aerosol droplets in different regions of the box model was measured with a Particle Counter Accumulator.

Parameter study
Given the lack of computational studies of the PIPAC technique, we investigated the influence of different parameters on aerosol droplet behavior. A baseline CFD model was generated using typical PIPAC parameter values (dp= 30 µm and Q= 0.5 mL s -1 ). [21] Subsequently, a parameter study was performed to vary the aerosol droplet diameter (dp [µm]), flow rate of nebulization (Q [mL s -1 ]), and liquid viscosity (µ [mPa s]). The simulations allowed to determine the effect of these parameters on the deposition of the aerosol droplets in different regions of the model.

Experimental validation of the CFD model
In a first step, we modelled filling of the box with CO2 gas to a pressure of 12 mmHg (1600 Pa), Figure  4. Next, we simulated the spatial distribution of aerosol droplets after nebulization for PIPAC and ePIPAC (Figure 5A and 5B). 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  5C 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. [27] As mentioned in section 2.3.6 (Equation 9), the Stk number , is a measure of the influence of the inertial effects during droplet transport. Thus, Stk as well as Rep 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 Equation (4), 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 3 rd power (see Equation 6). 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. [27] 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 Figure 6). The simulations showed that aerosol droplet deposition is significantly affected by their diameter. As shown in Figure 6, 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 7 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.

Effect of liquid flow rate
Fixed droplet diameter Figure 8A and 8B 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 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. [28] Also, low flow rates tend to decrease the nebulization angle (see section 3.5). 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 Figure 8A, 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 9 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.

Variable droplet diameter
The magnitude of inertial impaction depends on the droplet diameter and liquid flow rate. As proposed by Cheng et al., inertial impaction can be expressed as an impaction parameter, defined as dp 2 *Q, with dp=aerosol droplet diameter [µm] and Q=liquid flow rate [mL s -1 ] [29] . Table 1 shows impaction parameter values for six different flow rates of nebulization (Q=0.4, 0.5, 0.6, 0.7, 0.8 and 0.9 mL s -1 ). Figure 10 depicts a comparison of the simulated aerosol droplet deposition in the four regions of the CFD box model for different values of the impaction parameter (dp 2 *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 Figure 9, the clinical device generates a mean droplet diameter range between 25-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 viscosity
Dynamic viscosity measurements of Icodextrin resulted in values of 1.88 and 2.24 mPa s for concentrations of 4% and 7.5%, respectively. Keck et al. [30] performed a granulometric analysis for nebulized Icodextrin using Capnopen ® . Their results showed that the droplet size did not vary considerably with Icodextrin concentration: the median droplet diameter measured at two different Icodextrin concentrations (4% and 7.5%) was around 30 µm. Consequently, for the present CFD simulation, the droplet diameter was set at 30 µm with a flow rate of 0.5 mL s -1 . Figures 12 and 13 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.   Figure 14 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 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 15 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 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. In Figure 15, CFD modelling appears to give an underestimated value for proportion stained as compared to in vitro experiments. The lower estimates are probably explained by the fact that the Measurement Accumulator in COMSOL Multiphysics only considers completely black coverage of tissue, while dark gray was considered as well in vitro. 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.