Catheter placement selection for convection-enhanced delivery of therapeutic agents to brain tumors [version 1; peer review: awaiting peer review]

Background: Convection-enhanced delivery (CED) of therapeutic agents to brain tumors allows clinicians to bypass the blood-brain barrier (BBB) to infuse virus therapy, biological, or chemotherapy directly into a brain tumor through convection. However, the effectiveness of infusions via CED may depend on catheter placement. Methods: This study used diffusion maps from magnetic resonance imaging (MRI) of human brain tumors and computational fluid dynamics (CFD) simulations to assess therapy volume distribution percentages based on catheter placement locations. Results: The primary outcome showed differences in volume distribution based on the catheter placement location. Total tumor volume filled ranged from 144.40 mm 3 to 317.98 mm 3 . Percent filled of tumor volume ranged from 2.87% to 6.32%. Conclusions: The selection of the location for catheter placement using the region with the highest volume filled may provide optimal therapeutic effect. The researchers conclude that CFD may provide guidance for catheter placement in CED of therapeutic agents.


Introduction
When using intravenous therapy to treat brain tumors, the blood-brain barrier (BBB) appears to hinder the effectiveness of the therapy. The BBB consists of glial cells and pericytes, which enclose the endothelial layer of the blood vessels, and restricts fluid flow into the brain 1 . Glucose, electrolytes, vitamins, peptides, regulatory proteins, and ~2% of small-molecule drugs (<~400 Dalton (Da) with <~8 hydrogen bonds) can cross the BBB; however, large-molecule drugs and ~98% of smallmolecule drugs cannot adequately cross the BBB 2,3 . Therefore, clinicians may consider direct injection of therapeutic agents into the tumor instead of delivering the agent intravenously.
Within the brain, stand-alone diffusion is slow and therapeutic agents may migrate only 1 mm from the source in three days 5 . Stand-alone diffusion in simulations for our study caused therapy to migrate a maximum distance of 6.68 mm from the source. Convection (i.e. injection of therapeutic agents under pressure) can improve the distribution of the therapeutic agents, as opposed to diffusion, which is passive and occurs during molecule movement from an area of higher concentration to an area of lower concentration.
Bobo et al. 6 developed the CED method to increase therapeutic distribution in the brain. To bypass the BBB and deliver therapeutic agents directly to tumors, CED uses a catheter inserted into a cannula, which connects to a syringe pump. The pump introduces and maintains a pressure gradient, which activates mass fluid movement or convection. Using convection, Bobo et al. 6 observed at least a 100-fold increase in therapy distribution. However, this 100-fold increase did not necessarily translate to success in clinical trials.
Stine et al. 1 reviewed more than 15 trials and found that the outcomes using CED did not produce better results than standard treatment. For these trials, 10 clinical groups reported median survival, which ranged from 18.5 to 60 weeks 7-13 . These results suggest that either the therapy or the procedure needs improvement.
When modeling therapy distribution using CED, Støverud et al. 14 considered a region of the gray and white matter of the brain, not including tumor, along with porosity and patient-specific diffusivity and permeability parameters. They found that the therapy distribution followed the white matter region more than the gray matter region. Using tumor shapes that ranged from oblate to prolate, Sefidgar et al. 15 , found that larger tumors decreased therapeutic agent distribution in the interstitial fluid and that the prolate shape resulted in a wider range of distribution than other shapes. Zhan et al. 16 modeled distribution using a patient-specific tumor geometry and constant parameters for pressure, permeability, and diffusivity and found that the therapy type and infusion locations influenced the distribution. When infusing therapy from the tumor center, uniform distribution mostly occurred. However, infusion from peripheral locations seemed to limit spatial distribution closer to that location.
Recently Bhandari et al. 17 compared distribution of three different therapies within the brain tumor using constant diffusivity, but variable permeability and porosity based on tumor to contrast agent correlation. These researchers found that at the outset distribution appears greater in higher permeability regions, but then shifts and appears greater in higher porosity regions with more available space. Tumor cell density appeared to be an effectiveness indicator and noticeably varied by therapy and time. Conversely, Zhan et al. 18 observed that the permeability minimally affected the therapy distribution volume but did alter the distribution shape in space. CED appeared to cause the interstitial fluid pressure to increase near the infusion location.
The aim of this study is to assess therapy distribution in the tumor based on catheter placement location. Measurement of the volume distribution was a predictor of therapy effectiveness [19][20][21][22][23] .
Using CFD this study calculated therapy distribution percentages and total tumor volume filled with therapy. The study simulated transient therapy distribution in a patient-specific brain tumor using a pressure-based solver. The model analyzes the effect of convection and tumor properties such as geometry, diffusivity, permeability, porosity, and interstitial fluid pressure on therapy distribution within the tumor. T1-weighted imaging (T1W) and diffusion-weighted imaging (DWI) provided patient-specific geometry, diffusivity, and permeability information. Using patient-specific brain tumor characteristics may improve the accuracy of the volume distribution prediction. The researchers conducted sixteen simulations by dividing the tumor into four regions and within each region introduced therapy into four random locations.

Model development
Previously, the investigators, using a two-dimensional domain, detailed the image processing and simulation pre-processing stages, the order of convergence for the mesh, and the numerical scheme accuracy for this analysis 24 . In the model development, simulation results for a two-dimensional domain closely agreed with the analytical results and mesh spacing was 0.5% of the domain length 24 . Average mesh spacing for this phase of the study was approximately 0.64% of the domain length. An examination of results from the numerical schemes, which included the first order, second order, power law, quick scheme, and third order scheme suggested that higher order schemes compared favorably to the analytical results 24 .

Tumor geometry
The tumor geometry, which represents initial diagnosis, for this simulation resulted from an MRI brain exam obtained on a 3 Tesla Phillips Ingenia MRI machine (Phillips Ingenia scanner, Netherlands) at the University of Alabama at Birmingham Hospital. The use of patient-specific information requires a review by the Institutional Review Board (IRB) to determine human subjects research status. After reviewing the submitted application the IRB determined that this research is not human subjects research. The acquired sagittal T1W gradient-echo images ( Figure 1a) included slice thickness, repetition time, and echo time of 1.2 mm, 6.8 ms, and 3.3 ms, respectively. Axial diffusion tensor images resulted from the DWI of the same patient ( Figure 1b). The patient received the standard of care, which is fractionated radiation therapy and chemotherapy with the drug temozolomide at a dose of 0.1 mL/hour for 6 hours. Repetition time and echo time of the DWI were 7927.7 ms, and 70.0 ms, respectively. The b-value, in-plane resolution, and slice thickness were 800 s/mm 2 , 1.75 mm, and 2 mm, respectively, with 33 diffusion directions performed. Images used in this study are available as Underlying data 25 .
The stereolithography (stl) tumor geometry (see Figure 2) resulted from the Lesion GNB version 2 software 26 . Investigators utilized the Ansys version 19.1 ICEM meshing tool to generate the computational mesh from the stl 27 .

DWI to T1W transformation
The researchers used DSI Studio to transform the DWI onto the T1W 28,29 . The transformation, which produced the diffusion tensor and fiber tracking images (Figure 3), was necessary to ensure that the diffusion tensors aligned to the geometry spatial locations. MATLAB allowed further validation of the transformation with the spatial plotting of the geometry with the diffusion tensors overlaid (Figure 4) 30 . Validation can also be performed using GNU Octave.

Diffusion tensors
Variable tensors included Dxx, Dyy, Dzz (principal and main diagonal diffusion rates and directions) and Dxy, Dxz, Dyz (off-diagonal diffusion rates and directions) aligned with the x, represents concentration change with respect to time, u · ∇c represents convection, and ∇ · (D∇c) represents diffusion due to the concentration gradient.
A generation rate, defined as follows, was necessary to consider the effects of the infusion rate, source volume, and the density in the transport equation. S ϕ represents generation rate. IR represents infusion rate. S v represents source volume and ρ is density. The source represents the catheter placement location infused with one milliliter of an oncolytic herpes simplex virus for six hours 32 . Therefore, IR was 0.1 mL/hr. The minimum dose was 1×10 -10 mL/m 3 (1×10 -19 mL/mm 3 ).   Table 1.
four regions with four catheter placement locations per region. See Table 1 for input coordinates for each location. For each location, the radius was constant at 1.5 mm, which accounts for the radius of the catheter.

Simulations
Using Ansys, the researchers conducted 16 simulations of therapy distribution within the tumor 27 . Diffusivity and main diagonal   For the simulations, the solver type was pressure-based with absolute velocity formulation and transient solver time. The model was viscous and laminar. The solution method was a third-order spatial discretization. The stability condition for a numerical scheme determined an acceptable time step size that did not cause non-physical solutions. The following equation is useful for calculating the maximum time step size (M TSS ) and resolving the unsteadiness of the instance 27 .
where L scale (conservative) represents MIN(L vol , L ext ), U is maximum velocity at the domain boundary, L vol represents 3 V, V is the domain volume, L ext represents MAX(L x , L y , L z ), L x , L y , L z are domain extents (x,y,z direction), and ν is the kinematic viscosity. The M TSS and time steps were 10 seconds and 2,160 steps, respectively. Tumor porosity and average interstitial fluid pressure (IFP) of 0.6 and 266.65 Pa resulted from previous studies 33,34 .

Results
Infusion time for all simulations was 21,600 seconds (6 hours) with 0.6 mL of therapy infused per catheter. Simulations for four unique infusion locations per region resulted in therapy concentration ranges as depicted in Figure 6- Figure 9. The average velocity magnitude and pressure were 3.43×10 -7 mm/s and 443.66 Pa, respectively. Average distance from the source was 7.4 mm. The investigators ranked each location based on at least 1×10 -10 mL/m 3 (1×10 -19 mL/mm 3 ) distribution throughout the tumor (see Table 2).
In region one the location with maximum therapy distribution was A with 314.96 mm 3 or 6.26% and a maximum Euclidean distance from this location of 9.6 mm. Location B (-12.09 mm, 28.16 mm, 39.74 mm) represented the least therapy distribution at 213.83 mm 3 or 4.25%. In region two the location with maximum distribution was D with 317.98 mm 3 or 6.32% and least distribution for region three and was 209.81 mm 3 or 4.17%.
In region four the location with the maximum distribution was B with 224.90 mm 3 or 4.47% and a maximum distance from this location of 8.32 mm. Location D (-12.33 mm, 8.80 mm, 53.16 mm) in region four represented the lowest distribution percentage and was 2.87% or 144.40 mm 3 .
The tumor volume was 5031.35 mm 3 . The number of locations with volume filled greater than 250 mm 3 was six, while the locations less than 250 mm 3 was ten (see Figure 10).

Discussion
The primary outcome of this research showed differences in the therapy volume distribution based on the catheter placement location and suggested that location may influence the distribution and therapeutic value. Total volume filled by the therapy ranged from 144.40 mm 3 to 317.98 mm 3 . Percent filled of tumor volume ranged from 2.87% to 6.32%. The average velocity range of 2.45×10 -7 mm/s to 4.49×10 -7 mm/s caused the therapy to displace from the source by an average of 7.4 mm, which was reasonable based on an infusion time of six hours.
Previously, research groups defined therapeutic value using the ratio between specific therapy volume of distribution (V d ) and volume of infusion (V i ), which is the therapy plus carrier fluid 22,[35][36][37] . In vivo rodent and nonhuman primates' experiments allowed these researchers to determine V d by using image processing to measure the distribution of the specific therapy in the brain. Distribution for oncolytic viruses appeared to be unavailable because this distribution varies based on infusion location, tumor parameters, and therapy clearance 38 . The researchers in this study did not identify any in vivo V d results for the oncolytic virus. However, clinicians at the authors' institution are currently conducting a clinical trial to investigate immunotherapy in canines 39 . The trial is a regional collaboration designed to assess brain tumor therapies in humans and animals. Nevertheless, oncolytic viruses can reproduce in tumor cells destroying these cells without harming normal cells 32 .  Clinicians in the authors' institution select four catheter placement locations to maximize therapy volume distribution. If the clinicians select the four previously mentioned patient-specific locations from regions 1, 2, 3, and 4, the effectiveness of therapy may improve. The selection of locations within each region with the highest total volume filled or highest tumor percentage filled may provide the most optimal therapeutic value.
Researchers in this study analyzed 16 random locations, which provided a baseline mathematical prediction of the optimal catheter placement location. Clinicians may not currently use a mathematical model to select catheter placement locations, but instead select locations based on the avoidance of cell areas with visible signs of necrosis. This mathematical model may be an improvement over the current clinical method. In the next phase of this study the investigators will use a design optimization technique, which will allow the analysis of additional locations.

Conclusions
In this study, porosity and IFP resulted solely from previous studies; while diffusivity and permeability were mostly patientspecific. Patient-specific porosity and IFP may improve volume distribution accuracy. Although the authors did not consider retrograde infusion flow depending on tumor density and placement, it may improve the effectiveness of the treatment if future researchers consider the impact of retrograde flow.
The results presented suggest that computational fluid dynamic approaches using diffusivity and permeability parameters of actual patient data could greatly improve the treatment of adult and pediatric brain tumor patients by optimizing the placement of catheters in convection enhanced therapy. Using the specific anatomy of the patient, this novel method would optimize catheter placement to provide maximal tumor coverage of the therapeutic agent. This is the first report using diffusivity and permeability of real patient data and computational fluid dynamic modeling to guide catheter placement for convection enhanced delivery of a therapeutic agent. This predictive quantitative method to determine the ideal catheter placement location will assist clinicians in effectively treating brain tumors using CED.
This project contains T1-weighted and diffusion-weighted images used in the present study.
Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

Author contributions
Lisa H Antoine substantially contributed to the study design, data (T1-weighted images and diffusion-weighted images) acquisition, analysis and interpretation, drafting and critically revising the article, and the final approval of the published version. Roy P Koomullil substantially contributed to the study conception and design, data analysis and interpretation, critically revising the article, and the final approval of the published version. Timothy M Wick, Louis B Nabors, and Mark S Bolding substantially contributed to the study design, data analysis and interpretation, critically revising the article, and the final approval of the published version. Ahmed K Abdel Aal substantially contributed to data acquisition, analysis and interpretation, critically revising the article, and the final approval of the published version. and deformation during convection-enhanced drug delivery into brain tissue.