**1.Introduction**

Radioembolization is a treatment for patients with inoperable malignant liver tumors, one of the deadliest types of cancer worldwide [1]. The treatment is carried out using a microcatheter placed in the hepatic artery to infuse radiolabeled microspheres into the hepatic arterial bloodstream, which carries the microspheres to the tumoral bed, where they ge<sup>t</sup> lodged and irradiate tumoricidal doses of radiation [2]. The treatment outcome depends on the dose absorbed by tumors, meaning that it depends on the final distribution of microspheres in tumors. The absorbed tumoricidal dose to be achieved depends on the types of microspheres, which can be resin or glass microspheres loaded with isotope yttrium-90 or with isotope holmium-166 [3].

In the last decade, many computational fluid dynamics (CFD) studies have analyzed the microsphere distribution in hepatic arterial trees during radioembolization by analyzing the hepatic artery hemodynamics and microsphere transport, and assessed the influence of various parameters (e.g., injection velocity [4,5], microcatheter location [6–9], catheter type [10], etc.) on the microsphere distribution at the outlets of the trees. However, to date,

**Citation:** Lertxundi, U.; Aramburu, J.; Rodríguez-Fraile, M.; Sangro, B.; Antón, R. Computational Study of the Microsphere Concentration in Blood during Radioembolization. *Mathematics* **2022**, *10*, 4280. https:// doi.org/10.3390/math10224280

Academic Editor: Mauro Malvè

Received: 29 September 2022 Accepted: 14 November 2022 Published: 16 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

no study has focused on the analysis of the microsphere concentration in the blood (i.e., the number of microspheres per unit blood volume) as these microspheres exit the hepatic arteries under study. The analysis of this parameter could be included when developing an integrated software package to provide clinicians with practical recommendations about the optimal treatment parameters based on CFD simulations [11].

We hypothesize that the microsphere concentration in the blood could play a role in the distal penetration of microspheres and tumor microsphere coverage, and therefore in the final microsphere deposition in tumors. Indeed, recent publications have stressed the importance of glass microsphere density in tumors (i.e., the number of microspheres per unit tumor volume) [12,13], the number of microspheres injected based on the vascularity of tumors [14], and the method of administration [15] on the treatment outcome. Additionally, the microsphere volume fraction, which is directly related to the microsphere concentration in the blood, could play a role in the penetration of microspheres because the probability of clogging distal arterioles with vessel-to-microsphere diameter ratios below 3 depends on the microsphere volume fraction [16]. Such clogging may occur when microspheres would group together before reaching distal vessels with diameters below the microsphere diameter.

In addition to studying the microsphere distribution at the outlets of the arterial trees, CFD simulations allow for analyzing the microsphere concentration in the blood at such outlets. In order to explore this parameter, the aim of the present study was to analyze the influence of injection conditions on the microsphere concentration in the blood at the outlets of a patient-specific hepatic artery tree using CFD.

### **2. Materials and Methods**

### *2.1. Patient Data*

This study was based on a patient-specific case. To conduct the present study, the protocol 186/2018 was approved by the ethics committee of the University of Navarra and informed consent was signed by the patient. This patient was a 69-year-old male with multinodular hepatocellular carcinoma (HCC) involving liver segments S2, S3, S6, S7, and S8.

For CFD simulations, the geometries of the hepatic artery and microcatheter, and boundary conditions, were needed. The hepatic artery geometry and boundary conditions were extracted or derived from the information provided by MeVis (MeVis Medical Solutions AG, Bremen, Germany), that is, the three-dimensional hepatic artery, the healthy and tumor tissue volumes per liver segment, and information regarding the specific artery branches that feed each liver segment. In this case, the geometry of the hepatic artery starts at a 4.5 mm diameter proper hepatic artery (PHA) level and bifurcates until 43 outlets are obtained. Figure 1a shows the liver geometry, with the locations of tumors and the hepatic artery tree, and Figure 1b shows a model of the hepatic artery with the microcatheter, the location of which is indicated with an arrowhead. The hemodynamic conditions were derived based on volumetry analysis (volumes of normal and tumor tissue per segment) and perfusion-CT analysis (arterial perfusion of normal and tumor tissue), which were used to calculate the blood flow rate at the inlet of the PHA and the flow distribution in the 43 outlets [17]. This patient also participated in a previous study by the authors to validate the simulation methodology by comparing the microsphere distributions measured in vivo with the microsphere distributions calculated based on the results of CFD simulations that reproduced the actual treatments [18]. The segment-to-segment tissue volumes and arterial blood flow rates are tabulated in Table 1. Even though, previously, three microsphere injections were administered to the patient, in this study, only one injection was simulated, with the microcatheter placed in the PHA, 20 mm away from the bifurcation. The microcatheter was modeled as a standard end-hole microcatheter with an inner diameter of 0.65 mm and outer diameter of 0.9 mm.

**Figure 1.** (**a**) MeVis study showing the liver, hepatic arterial tree, and tumor nodules. (**b**) Hepatic artery tree model with the microcatheter tip indicated by an arrowhead. PHA: proper hepatic artery; LHA: left hepatic artery; RHA: right hepatic artery.


**Table 1.** Patient's liver mass volumes and blood flow rates per segment.

### *2.2. Simulation of Radioembolization*

The geometry of the hepatic artery was imported into the software package SpaceClaim (Ansys Inc., Canonsburg, PA, USA), where the microcatheter geometry was also added. The geometries were discretized in Fluent Meshing 2021R1 (Ansys Inc.) using poly-hexcore elements, yielding a sufficiently fine mesh of 2.5 million elements.

Regarding the governing equations and boundary conditions, the model used by Aramburu et al. [6] was used: blood flow was simulated using the conservation of mass and the conservation of linear momentum in a laminar regime for an incompressible (1050 kg/m3) and non-Newtonian fluid with the viscosity modeled using a modified Quemada model. Microspheres were modeled as 32 μm, 1600 kg/m<sup>3</sup> spheres and their dynamics were modeled with Newton's second law of motion. Four forces were considered to be acting on the microspheres: virtual mass force, gravitational force, pressure gradient force, and drag force. The interaction of blood flow with microspheres was modeled as bidirectional, but the interaction between microspheres was neglected. As for the boundary conditions, a periodic pulsatile fully-developed velocity profile with a period of one second representing one cardiac cycle was prescribed at the inlet, flow fractions at the outlets, and the no-slip condition at the walls. Microspheres were injected through the inlet of the microcatheter and elastic collisions were assumed at the walls.

The equations of the model were solved numerically using Fluent 2021R1 (Ansys Inc.). The SIMPLE scheme (i.e., Semi-Implicit Method for Pressure Linked Equations) was used for the pressure and velocity coupling, the least squares cell-based algorithm for the computation of gradients, and second-order schemes to interpolate pressure and momentum. The convergence criterion was to attain scaled residuals of 10−5. The time step was fixed at a value of 2 milliseconds with a maximum of 80 iterations per time step. A total of seven cardiac cycles were simulated: the first cycle, i.e., one second, was for flow convergence (from *t* = −1 s to *t* = 0 s), during the second cycle the microspheres were injected (from *t* = 0

s to *t* = 1 s), and the remaining five cycles ensured that the majority of the injected microspheres exited the domain (from *t* = 1 s to *t* = 6 s). Simulation results were analyzed from *t* = 0 s onward. In this study, depending on the simulation, the injection of microspheres sometimes lasted more than one cycle, in which cases the remaining cycles for microsphere exiting were fewer than five cycles.

### *2.3. Study Design and Postprocessing*

The objective of this study was to analyze the influence of injection conditions on the microsphere concentration in the blood at the outlets. The number of injected microspheres can be calculated using Equation (1):

*N* = *<sup>C</sup>*vial*U*inj*AcT*, (1)

where *N* (microspheres, hereafter MS) is the total number of microspheres injected in the simulation, *C*vial (MS/mL) is the number of microspheres per unit milliliter in the vial, *<sup>U</sup>*inj (cm/s) is the injection velocity, *Ac* (cm2) is the cross-sectional area of the microcatheter, and *T* is the injection span (s).

In this study, three simulations were run. In these simulations, *N* was kept constant, meaning that the same number of microspheres was injected, and *<sup>U</sup>*inj was also kept constant, so *C*vial and *T* were modified, as shown in Table 2.

**Table 2.** Characteristics of the simulations of the study. In all simulations, the number of injected microspheres and injection velocity were kept constant; in Simulation #2, the injection span and the concentration of microspheres in the vial were twice and half of the injection span and the concentration of microspheres in the vial in Simulation #1, respectively; in Simulation #3, the injection span and the concentration of microspheres in the vial were three times and one-third of the injection span and the concentration of microspheres in the vial in Simulation #1, respectively.


Regarding the postprocessing of simulation results, the overall outlet-to-outlet and segment-to-segment microsphere distributions were extracted, as well as the concentration of microspheres in the blood in the PHA and at the outlets over time. To calculate the concentration of microspheres in the blood in the PHA, Equation (2) was used:

$$\mathcal{L}\_{\text{PHA}} = \frac{N\_{\text{PHA}}}{V\_{\text{PHA}}},\tag{2}$$

where *C*PHA (MS/mL) is the concentration of microspheres in the blood in the PHA calculated for a given span, *N*PHA (MS) is the number of microspheres that were infused into the PHA during the same span, and *V*PHA (mL) is the volume of blood and the volume of injection fluid that flowed through the PHA during the same span. To calculate the concentration of microspheres in blood at the outlets Equation (3) was used:

$$C\_i = \frac{N\_i}{V\_{\text{blood},i}},\tag{3}$$

where *Ci* (MS/mL) is the concentration of microspheres in the blood for outlet *i* (with *i* from 1 to 43), calculated for a given span, *Ni* (MS) is the number of microspheres that exited through outlet *i* during the same span, and *<sup>V</sup>*blood,*i* (mL) is the volume of blood that flowed through outlet *i* during the same span. In this study, a span of 20 milliseconds was selected.
