Magnetic resonance (MR) images of a healthy adult subject were used to create realistic FE model of the brain in 3D. The MR images were acquired from our previous study (Angeli and Stylianopoulos, 2016) and were reused in this work. Specifically, for the morphological imaging of the brain a T1-weighted, three-dimensional, fast field echo pulse sequence was acquired with an echo and repetition time of 3.2 and 7.1 ms respectively, while an isotropic voxel size of 1 mm was used to cover the entire brain. The commercial software ScanIP from Simpleware (version 6.0; Synopsys, Mountain View, United States) was employed for the three-dimensional reconstruction of the brain geometry. Specifically, two masks were first generated from the MR images using Simpleware’s “threshold” operation, which selects each pixel according to its brightness. The darkest areas comprise the mask of the gray matter, whereas the brightest regions comprise the mask of the white matter in the brain. These masks were the two different domains of the resulting 3D brain geometry. Then, the “island removal” and “cavity fill” operations were used to eliminate small unconnected parts of the masks and fill any gaps of the model, respectively. Additionally, smoothing was performed using Simpleware’s “Gaussian smoothing” operation. The 3D brain geometry was eventually created and exported in a COMSOL-compatible geometry file. The 3D brain geometry was imported to the commercial FE software COMSOL Multiphysics (version 5.5; COMSOL Inc., Burlington, MA, United States). Subsequently, a sphere with radius 6 mm and a cylinder with radius 1.5 mm were formed inside the brain geometry to represent the tumor and catheter domains, respectively. The sphere (brain tumor) was located in the white matter inside the parietal lobe of the left cerebral hemisphere. Hence, the 3D brain geometry consisted of four domains, i.e., gray matter, white matter, tumor, and catheter. Next, using COMSOL Multiphysics software, the FE mesh was created though the “free tetrahedral” and “boundary layer” operations to form an optimal mesh with boundary layers at the interfaces between the geometry domains. Figure 1 illustrates the individualized 3D FE brain model used for the CED simulations.

The brain tissue was modeled as a porous medium (Angeli and Stylianopoulos, 2016) and the drug distribution and pharmacodynamics model was based on the conservation equations of mass transport and momentum balance. The drugs were released directly into the tumor interstitial space through the catheter and then they could travel inside the tumor by convection and diffusion, exit the tumor region towards the surrounding healthy tissue and be cleared by tumor and healthy vessels by convective and diffusive mechanisms.

Interstitial flow within the porous brain tissue obeys to the continuity equation and the extended Darcy’s law (Nield and Bejan, 1999). The continuity equation of tumor fluid phase is given by: ∇ ⋅ u = Q m , (1) where u is interstitial fluid velocity (IFV) and Q _m is the fluid flux that is exchanged between the tissue and the vascular system, and it is mathematically modelled using Starling’s approximation (Vavourakis et al., 2017; Lambride et al., 2020): Q m = L p S v ( p v − p i ) − L p l S v l ( p i − p l ) , (2) where L _p , S _v, and p _v are the hydraulic conductivity of the vessel wall, the vascular density, and the vascular pressure of the blood vessels respectively, whereas L _pl , S _vl , and p _l are the corresponding quantities for the lymphatic vessels, and p _i is the interstitial fluid pressure (IFP). The first term in Eq. 2 refers to the fluid flux entering the tumor or the surrounding normal tissue from the blood vessels, while the second term refers to the fluid flux exiting through the lymphatic system. It is important to note that the second term was set to zero in the tumor tissue due to the dysfunctional lymphatic system (Stylianopoulos et al., 2018).

We extended Darcy’s law due to the presence of infusate flow from the CED catheter, to describe the fluid motion in porous tissues of the brain (Nield and Bejan, 1999; Linninger et al., 2008). The Brinkman equation extends Darcy’s law to describe the dissipation of the kinetic energy by viscous shear, similar to the Navier-Stokes equation. Thus, the momentum balance is given by: ρ ∂ u ∂ t = ∇ ⋅ [ μ i ε p ( ∇ u+u ∇ ) − ( 2 3 μ ε p ∇ ⋅ u + P i ) I ] − ( 1 k + Q m ε p 2 ) u , (3) where ε _p is the tissue porosity, k is the hydraulic conductivity of the interstitial space, μ _i and ρ refer to the viscosity and the density of interstitial fluid, respectively.

Regarding the drug concentration in the porous brain tissue, the convection-diffusion mass transport equation was applied, incorporating the exchange of drug between tissue and the blood vessels with the term Q (Mpekris et al., 2017): ∂ C i ∂ t + ∇ ⋅ ( C i u − D i ∇ C i ) = Q − k d C i − k l C i , (4) where C _i refers to the relative concentration of drug in the interstitial space, which is a dimensionless quantity, i.e., the ratio of the local drug concentration to the value of the concentration entering the catheter, and D _i is the diffusion coefficient of the drug in the interstitial space. Additionally, k _d is the drug degradation rate constant, while K l is the lymphatic drainage constant, applicable only in the healthy tissue region (Linninger et al., 2008; Stylianopoulos et al., 2018).

The total net transvascular drug flux between the brain tissue and the vessels can be written neglecting the oncotic pressure difference across the wall (Mpekris et al., 2017; Stylianopoulos et al., 2018): Q = ( C v − C i ) P S v + ( 1 − σ f ) ( p v − p i ) L p S v C v , (5) where p is the vascular permeability of the drug through the pores of the vessel wall, C _v is the vascular concentration of the drug, and σ _f is the reflection coefficient. In CED, the therapeutic agents are injected directly inside the brain tumor tissue through a catheter; thereby in the tumor domain the vascular drug concentration (C _v ) is considered negligible compared to the interstitial space drug concentration (C _i ).

The therapeutic agent was considered to be of spherical shape and without any surface charge and the vessel wall openings were modeled as cylindrical pores (Deen, 1987). To investigate the direct effect of the vascular wall pore size and drug size on the drug concentration distribution, the theory for hindered transport of rigid solutes through liquid filled porous was introduced. Using this theory, the hydraulic conductivity of vascular walls, L _p , the drug vascular permeability, p, and the reflection coefficient, σ _f , were explicitly estimated based on the ratio of particle size to the vessel wall pore size. The hydraulic conductivity of vessel walls is calculated via Eq. 6 and the vascular permeability of the drug is defined by Eq. 7: L p = γ r 0 2 8 μ L v w , (6) P = γ H D L v w , (7) where γ is the fraction of vessel all surface area occupied by pores, r _o is the pore radius, μ is the viscosity of plasma at 310 K, L _vw is the thickness of the vessel wall, and H corresponds to the hydrodynamic coefficient for neutral spheres in cylindrical pores. D is the diffusion coefficient of a particle with radius r _s , given by the Stokes-Einstein relationship, D = K b T 6 π μ r s , where K _b is the Boltzmann constant, and T is the temperature of the solution (Deen, 1987). Definition of the reflection coefficient, σ _f , and the hydrodynamic coefficient, H, is provided in Supplementary Table S2.

For a realistic 3D finite element (FE) brain model, data from experimental studies were used to determine the physiological and biomechanical properties. As in pertinent studies, the diameters of the vessel wall pores were set to 50, 100, and 150 nm (Hobbs et al., 1998; Sarin et al., 2009; Chauhan et al., 2012), whereas the hydraulic conductivity of the brain tumor interstitial space, k, has been reported to vary significantly, as a result it was taken equal to 2 × 10⁻¹⁴, 2 × 10⁻¹³, and 2 × 10⁻¹² m²/(Pa s) (Netti et al., 2000; Smith and Humphrey, 2007). Optimal size of nanoparticles for cancer treatment ranges from 20 to 70 nm in diameter to easily diffuse into tumor ECM (Koo et al., 2006; Vavourakis et al., 2018; Wijeratne and Vavourakis, 2019). Recent studies have showed that particles of sizes 40–50 nm are able to effectively bind and induce receptor-mediated endocytic processes (Jain and Stylianopoulos, 2010). Additionally, it has been shown that particles larger than 60 nm in diameter are less effective diffusing through the extracellular matrix space (Jain and Stylianopoulos, 2010). With this in mind, we decided to focus on the drug-size range of 1–60 nm in diameter referring to conventional chemotherapeutics (<2 nm), antibodies and liposomes. For a relevant drug-size comparison summary for brain tumors, the reader should refer to Table 1 from (Koo et al., 2006). Thus, the effective diameter of the therapeutic agents varied from 1 nm for small-size molecules, to 20 and 60 nm for liposomes to investigate the delivery of a wide range of drug sizes. All values of the model parameters are summarized in Supplementary Table S3. It is noteworthy that for patient-specific studies, the parameters in Eq. 5 can be estimated for a small drug by direct measurement in a DCE imaging protocol to accurately define the drug capillary loss rate (Brady et al., 2020).

For each geometry domain, the conservation equations of mass and momentum were discretized and solved numerically coupled using the commercial FE software COMSOL Multiphysics 5.5 (COMSOL Inc., Burlington, MA, United States). The Brinkman equation module of the FE software was used, providing the conservation of mass and momentum that can fully describe the fluid dynamics within each geometry domains. The Brinkman equation module computes both the velocity filed and pressure, which are the dependent variables. Additionally, both the velocity (vector) field and pressure (scalar) were discretized using linear Lagrange basis functions. Likewise, the drug concentration (scalar) was discretized using linear Lagrange basis functions. COMSOL project files for two representative simulation cases of CED can be freely accessed on Figshare: https://figshare.com/projects/Convection-Enhanced_Delivery_in_silico_study_for_personalized_brain_cancer_treatment/135923.

Following data from relevant clinical CED studies (Barua et al., 2013; Shi and Sanche, 2019; Stine and Munson, 2019; Tosi and Souweidane, 2020), the infusion flow rate, Q _f , ranged from 0.025 to 0.75 ml/h and the infusion volume varied between 0.25 and 185 ml. A flow rate of 0.5 ml/h and an infusion volume of 3 ml were selected for all simulations. Therefore, at the interface between the catheter and tumor the normal inlet velocity (i.e., U _o =Q _f /A, where A the cross section of the catheter) was taken equal to 1.99 × 10⁻⁵ m/s and infusion lasted for 6 h. Also, at the interface of the catheter and tumor tissue, the relative drug concertation was set to unity for the period of the infusion and after completion of infusion, a zero-flux boundary condition was applied (i.e., n ⋅ ( D ∇ C i − u C i ) = 0 , where n corresponds to the outward unit normal vector).The normal stresses on the outer brain surfaces were equal to zero (i.e., n σ = 0). At the catheter surfaces, a no-slip boundary condition was applied for fluid velocity (i.e., u = 0) and additionally, a zero-flux boundary condition was set for the transport of the drugs. The latter boundary condition was also set at the outer brain surfaces (Supplementary Figure S1). Regarding initial values, both the fluid velocity and pressure variables were set to zero for all geometry domains. Likewise, the initial value of the drug concentration was zero everywhere in the brain model.

We first set the hydraulic conductivity of the tumor interstitial space equal to 2 × 10⁻¹³ m²/(Pa s) and the diameter of vessel wall pores to 100 nm (baseline tumor microenvironment conditions). Figure 2 illustrates a sagittal view of the center of the tumor tissue, presenting the spatial distribution of different drug sizes in three snapshots, i.e., 6, 12, and 24 h, respectively after CED injection commenced. All spatial drug distributions are symmetric in the vertical axis. In case of the 1 nm drug size diameter, the highest drug concentration during infusion is located near the infusion site. Small molecule therapeutics owing to their small size have a high diffusion coefficient and thus, they travel fast away from the infusion site and can be cleared easily by the pores of the blood vessels. Hence, a significant drug amount is lost because the drug can easily pass into the blood vessels due to concentration gradients. These observations for small drug sizes were also reported by (Brady et al., 2020). Brady et al. (2020) mentioned that the capillary loss rate plays a key role in the drug transport within the tumor tissue, especially for small molecules in which the loss rate is significantly high. As a result, the concentration of the drug in the tumor diminishes soon after the infusion stops (Figure 2A). When the infusion stops, the average relative concentration of the therapeutic agent reaches 0.25 in the tumor tissue, and then decreases sharply to zero (Figure 3A).

For the 20 nm therapeutic agent, vascular permeability becomes 98% lower compared to that of the 1 nm drug size, and thereby, the intravasation of the drug is hindered significantly. This type of nanoparticle is still able to diffuse relatively fast and penetrate the tumor tissue, covering almost the entire tumor region during CED administration (Figure 2B). The maximum average relative concentration of drug in the tumor tissue is 0.78 at the end of the infusion period at 6 h but drug remains in the tumor region for over 24 h (Figures 2B, 3A). Similar to the 20 nm drug size, the 60 nm drug is distributed in a large volume of tumor tissue during infusion (Figure 2C). The maximum average relative concentration in the tumor tissue is equal to 0.81 at 6 h. Additionally, 6 h after initial injection, the drug is dispersed in the same way as the 20 nm drug but the average relative concentration of drug in the tumor is higher during the whole period (Figure 3A).

The hydraulic conductivity of the tumor describes the resistance to interstitial fluid flow through the pores of the interstitial space of the tissue and thus, it is directly related to IFV and the convective transport of drugs. We repeated simulations to investigate the effect of hydraulic conductivity of the tumor interstitial space on drug concentration and spatial distribution. To achieve this, the hydraulic conductivity was decreased from 2 × 10⁻¹³ to 2 × 10⁻¹⁴ m²/(Pa s), whereas the pore diameter of tumor vessel walls was kept to the baseline value of 100 nm. This modification affects insignificantly the spatial distribution of the 1 nm drug (Figures 4A, 3B). In the case of the 20 and 60 nm drugs, higher drug concentrations are observed in the tumor center after infusion (Figure 4A). This can be justified by the IFP and IFV profiles (Supplementary Figure S2). As the hydraulic conductivity of the tumor interstitial space decreases, the velocity within the tumor decreases, while the fluid pressure increases in the tumor center. Therefore, drug transport through convection is reduced in the tumor center, rendering diffusion the dominant transport mechanism, which is inversely proportional to the size of the drug. The convection contribution increases near the tumor boundary, where the drug concentration is low. The maximum drug concentration is located both in the center of the tumor tissue and near the catheter outlet and the drug spreads sufficiently within the tumor tissue. For 20 and 60 nm drug diameter, the average relative concentration values after 12 h are 0.32 and 0.39, respectively, ∼78% greater with respect to the baseline conditions (Figure 3B).

Subsequently, we increased the value of the tumor tissue hydraulic conductivity to 2 × 10⁻¹² m²/(Pa s), an order of magnitude higher than baseline value. By increasing the hydraulic conductivity of the tumor tissue, the IFP decreases and the IFV rises within the tumor tissue. Thus, convective drug transport in the tumor interstitium is enhanced from the center towards the tumor periphery, allowing the drug to escape easily from tumor tissue. Small molecule drugs with diameters up to 1 nm behave similar to previous simulations as they still diffuse fast, and their transport does not depend on pressure gradients. However, after CED administration, the concentration of all drug sizes is at low levels inside the tumor (Figure 4B). The distributions of nanoparticles with 20 and 60 nm diameters are identical, and the average concentrations decrease sharply as a function of time (Figure 3C). Specifically, after 12 h, the average relative concentration values in the tumor for 20 and 60 nm drugs are 0.048 and 0.056 respectively, ∼74% lower with respect to the baseline conditions.

Next, we investigated the effect of the pore size of the tumor vessel walls on the distribution of drugs administered through CED. A series of simulations was performed varying the pore diameters of tumor vessel walls and calculated the spatiotemporal distribution of the drugs. The drug diameter range remained the same with the previous simulations and the hydraulic conductivity of tumor tissue was set to 2 × 10⁻¹³ m²/(Pa s). By changing the pore size of vessel walls from 50 to 150 nm, the hydraulic conductivity of the vessels increases considerably, according to Eq. 6. As a result, the fluid pressure increases uniformly inside the tumor and drops steeply at the tumor margin, whereas the velocity magnitude increases at the periphery following the IFP gradients (Supplementary Figure S3), which is a hallmark of tumor patho-physiology (Baxter and Jain, 1989). Taken together, this means that increasing the pore diameter of the vessel walls reduces drug transport through convection. Therefore, the main transport mechanism is through diffusion. Comparing the drug distribution with the baseline tumor microenvironment conditions (Figure 2), the modification of pore dimeters does not influence the spatial distribution of small drugs up to 1 nm diameter (Figure 5). However, the distribution of nanoparticles of 20 and 60 nm diameter differ significantly with respect to the tumor vessel pore size. In the case of 50 nm diameter endothelial wall pores, after infusion, the drug disperses homogeneously within the tumor tissue due to enchanted convection (Figure 5A). Furthermore, it is observed higher local drug concentration for 150 nm pore diameter (Figure 5B). After 12 h, as shown in Figure 6, the average relative concentration for 60 nm drug size is 0.25% ∼14% greater with respect to the baseline conditions.

Supplementary plots for each drug size presenting the average relative drug concentration as a function of time for different values of the hydraulic conductivity of the tissue interstitial space and for different diameters of vessel wall pores are presented in Supplementary Figures S4–S6. These comparative plots show separately the effect of tumor microenvironment properties on the average relative concentration for each drug.

To investigate the heterogeneity of the drug distribution, the drug concentration as a function of distance from the tumor center was calculated for four different directions in the plane (Figure 7A) and the mean value as well as the standard error of the results along these four directions are presented (Figures 7B–E). The standard error bars represent the heterogeneity of the drug distribution, giving a quantification of how close the results are to the average value, i.e., variation of average value. The smaller the error bars, the more homogeneous the distribution of the drug is. The heterogeneity of the drug distribution was assessed during and after CED administration, taking snapshots at 5 and 9 h. The drugs with a diameter of 20 and 60 nm were selected for the heterogeneity measurements because the concentration of the 1 nm drug decreases sharply to zero after infusion. Also, the baseline tumor microenvironment conditions, i.e., 100 nm pore diameter of tumor vessel walls and a 2 × 10⁻¹³ m²/(Pa s) hydraulic conductivity of tumor interstitial space were assumed.

According to the error bars, there is a homogeneous spatial distribution of the drug near the tumor center for both drug diameters during infusion. The heterogeneity exists radially outwards where the magnitude of the error bars increases (Figures 7B,D). This happens because lines 1 and 2 are closer to the infusion site and the drug concentration values are higher along their lines compared to lines 3 and 4 (Figure 7A). After infusion, the average concentrations for both drug sizes gradually decrease as a function of the distance and the spatial distribution of the drug becomes homogeneous (Figures 7C,E), as it is also verified from Figures 2B,C. Furthermore, we repeated simulations varying the tumor hydraulic conductivity (Supplementary Figures S7, S8) with the results following the same patterns as in Figure 7.

Finally, we investigated the misplacement of the CED catheter. In contrast to the previous simulations where the catheter injects the drug directly into the tumor tissue, the following series of simulations were performed to examine the case that the catheter was taken to be outside of the tumor tissue. To investigate the effect of the catheter placement on drug delivery, the CED catheter was placed 1 mm outside the tumor, as illustrated in Figure 8. During infusion, drug delivery into the tumor is hindered by the pressure gradients at the tumor/heathy tissue interface that drives fluid flow from the tumors towards the healthy tissue. Therefore, high concentrations are found in the healthy tissue. In the case of drug with 1 nm, a tumor volume percentage for relative concentration above 0.2 is approximately 7.5% at the end of the infusion period. For larger drug sizes, after infusion, drug concentration remains high in the tumor tissue. Specifically, after 12 h, the upper part of tumor tissue has 0.8 relative drug concentration and the tumor volume percentage for relative concentration above 0.2 is 12.5% for the 60 nm therapeutic agent.