Convection-Enhanced Delivery of Antiangiogenic Drugs and Liposomal Cytotoxic Drugs to Heterogeneous Brain Tumor for Combination Therapy

Simple Summary Although developed anticancer drugs have shown desirable effects in preclinical trials, the clinical efficacy of chemotherapy against brain cancer remains disappointing. One of the important obstacles is the highly heterogeneous environment in tumors. This study aims to evaluate the performance of an emerging treatment using antiangiogenic and cytotoxic drugs. Our mathematical modelling confirms the advantage of this combination therapy in homogenizing the intratumoral environment for better drug delivery outcomes. In addition, the effects of local microvasculature and cell density on this therapy are also discussed. The results would contribute to the development of more effective treatments for brain cancer. Abstract Although convection-enhanced delivery can successfully bypass the blood-brain barrier, its clinical performance remains disappointing. This is primarily attributed to the heterogeneous intratumoral environment, particularly the tumor microvasculature. This study investigates the combined convection-enhanced delivery of antiangiogenic drugs and liposomal cytotoxic drugs in a heterogeneous brain tumor environment using a transport-based mathematical model. The patient-specific 3D brain tumor geometry and the tumor’s heterogeneous tissue properties, including microvascular density, porosity and cell density, are extracted from dynamic contrast-enhanced magnetic resonance imaging data. Results show that antiangiogenic drugs can effectively reduce the tumor microvascular density. This change in tissue structure would inhibit the fluid loss from the blood to prevent drug concentration from dilution, and also reduce the drug loss by blood drainage. The comparisons between different dosing regimens demonstrate that the co-infusion of liposomal cytotoxic drugs and antiangiogenic drugs has the advantages of homogenizing drug distribution, increasing drug accumulation, and enlarging the volume where tumor cells can be effectively killed. The delivery outcomes are susceptible to the location of the infusion site. This combination treatment can be improved by infusing drugs at higher microvascular density sites. In contrast, infusion at a site with high cell density would lower the treatment effectiveness of the whole brain tumor. Results obtained from this study can deepen the understanding of this combination therapy and provide a reference for treatment design and optimization that can further improve survival and patient quality of life.


Introduction
Brain tumors have posed a unique health concern worldwide. Glioblastoma is the most aggressive primary brain cancer that remains incurable. The median survival is where t is time. ρ is the interstitial fluid density, and µ is its dynamic viscosity. The interstitial fluid pressure (IFP) and velocity (IFV) are represented by p ISF and v, respectively. κ is the tissue permeability that reflects the resistance of tissue to the interstitial fluid flow. F BL is the flux of fluid gain from the blood (BL), driven by the effective transvascular pressure gradient. It can be calculated by Starling's law, as where L BL is the hydraulic conductivity of the capillary wall. S BL /V T = ϕS BL,0 /V T is the area of capillary surface per unit tissue volume, with the subscript 0 denoting the initial condition. ϕ is the scaling factor that reflects the antiangiogenic effect. p BL is the blood pressure. σ BL is the osmotic reflection coefficient for proteins in the blood. π BL and π ISF are the osmotic pressure of blood and interstitial fluid, respectively. A brain, including the embedded tumor, consists of the microvasculature network and tissue; the latter can be further divided into the extracellular space (ECS), cell membrane (CM) and intracellular space (ICS). The transport of infused drugs between these compartments is schematically shown in Figure 1, where the letters LP, FD and BD stand for the liposome-encapsulated drugs, released free drugs and the drugs that bind with proteins, respectively. By assuming the liposomes are stable before entering a tumor, the concentration of liposome-encapsulated drugs (C LP,ECS ) can be described by where ECS is the volume fraction of tissue extracellular space. D LP,ECS is the diffusivity of liposomes in tissue ECS. k rel is the drug release rate from liposomes. k LP,b = P TV,LP S BL V T is the loss rate of liposomes to the blood, where P TV,LP is the liposome transvascular permeability.
The concentration of free drugs (C FD,T ) and bound drugs (C BD,T ) in the tissues can be accounted by where BL , CM and ICS are the volume fractions of plasma, cell membrane and intracellular space, respectively. C BD,CM is zero as no drugs would be eliminated or bound on the cell membrane [22]. The accumulation of free drugs in tissue is determined by multiple factors, including the dynamic release from liposomes, transport by convection and diffusion, blood drainage, drug elimination due to physical degradation and bioreaction, binding with proteins and cell uptake. Based on the law of conservation of mass, the availability of free drugs in tissue is governed by Drug transport in brain tumor after infusion. Liposomal cytotoxic drugs and antiangiogenic drugs are infused into the tumor ECS simultaneously. The released free cytotoxic drugs can cross the cell member to enter the cell interior; whereas liposome cell uptake strongly depends on the formulation, morphology and surrounding environment. As a result, some types of liposomes can enter cells efficiently, while others cannot cross the cell membrane. Therefore, this process is marked in a dashed line. Since there is a lack of a general model to govern this highly liposomespecific process due to its complexity, liposome cell uptake was ignored in several past studies [15,37]. AA-antiangiogenic drugs, LP-liposome-encapsulated cytotoxic drugs, FD-free cytotoxic drugs, and BD-cytotoxic drugs that bind with proteins.
The accumulation of free drugs in tissue is determined by multiple factors, including the dynamic release from liposomes, transport by convection and diffusion, blood drainage, drug elimination due to physical degradation and bioreaction, binding with proteins and cell uptake. Based on the law of conservation of mass, the availability of free drugs in tissue is governed by where FD,ECS is the diffusivity of free drugs in tissue ECS.  [39,40]. Therefore, Equation (6) can be rewritten as where = ECS (1 + BC ) + ICS IE (1 + BC ) + CM CE is a parameter determined by the properties of the drug and tissues.
The concentration of antiangiogenic drugs in the tissue ECS depends on the convective and diffusive transport, and elimination, as Drug transport in brain tumor after infusion. Liposomal cytotoxic drugs and antiangiogenic drugs are infused into the tumor ECS simultaneously. The released free cytotoxic drugs can cross the cell member to enter the cell interior; whereas liposome cell uptake strongly depends on the formulation, morphology and surrounding environment. As a result, some types of liposomes can enter cells efficiently, while others cannot cross the cell membrane. Therefore, this process is marked in a dashed line. Since there is a lack of a general model to govern this highly liposomespecific process due to its complexity, liposome cell uptake was ignored in several past studies [15,37]. AA-antiangiogenic drugs, LP-liposome-encapsulated cytotoxic drugs, FD-free cytotoxic drugs, and BD-cytotoxic drugs that bind with proteins.
where D FD,ECS is the diffusivity of free drugs in tissue ECS. k FD,b is the loss rate of free drugs to the blood, calculated by k FD,b = P TV,FD S BL V T . P TV,FD is transvascular permeability of free drugs. k FD,e is the elimination rate of free drugs, integrating the contributions of physical degradation and bioreaction. Two assumptions are further introduced: (a) The concentration of free drugs and bound drugs are linearly correlated in the tissue ECS and ICS (C BD,ECS = K BC C FD,ECS ; C BD,ICS = K BC C FD,ICS ) [23,38]; (b) Dynamic equilibrium is reached for the concentrations of free drugs between the tissue compartments (C FD,ICS = H IE C FD,ECS ; C FD,CM = H CE C FD,ECS ) [39,40]. Therefore, Equation (6) can be rewritten as where ω = ECS (1 + K BC ) + ICS H IE (1 + K BC ) + CM H CE is a parameter determined by the properties of the drug and tissues. The concentration of antiangiogenic drugs in the tissue ECS depends on the convective and diffusive transport, and elimination, as where D AA,ECS is the diffusivity of antiangiogenic drugs in the tissue ECS. k AA,e is its elimination rate. The antiangiogenic effect can be described by in which ϕ is the scaling factor for microvasculature density. k a is the antiangiogenic rate. α, β and γ are applied to consider the natural angiogenesis of the tissue.

MR Imaging and Processing Protocol
MR scan of two brain tumor patients was performed on a 3.0-Tesla Ingenia MRI scanner (Philips Healthcare, Amsterdam, The Netherlands) at Fortis Memorial Research Institute, Gurugram, India. MR imaging protocol included a standard protocol for brain tumor patients along with DCE MRI or T1-Perfusion MRI. DCE MR images of the brain were acquired using the 3D fast field echo (T1-FFE) sequence before, during and after contrast injection. In the current study, a systemic infusion of the contrast agent Dotarem (Gadoterate meglumine, Guerbet, France) with a dose of 0.1 mmol/kg body weight at an infusion rate of 3.0 mL/s was used. The key imaging parameters are summarised in Table 1. All procedures performed in this study involving human participants were in accordance with the ethical standards of the Institutional Review Board (IRB No. 2020-001- [19][20][21][22][23][24][25][26][27][28] and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. The demographic characteristics for the two exemplary patients are given in Table 2, where a pre-2016 classification is applied.  Figure 2. The brain tumor is segmented from its holding tissue based on the local signal intensity on each image slice. These segmentation results are then stacked along the MR scan direction to reconstruct the brain tumor in 3D. The tumor equivalent radius and volume are measured as 16 mm and 5461 mm 3 , respectively. Given this study is focused on the drug transport in the brain tumor rather than the whole brain, a rectangular volume [41,42] of interest with the dimension of 58 × 97 × 72 mm 3 is generated to fully enclose the brain tumor.
The final computational mesh consists of 82,944 uniform structured elements, obtained after the mesh-independence test [42]. The dimension of each element is 0.9 × 0.9 × 6.0 mm. Drugs are infused into the brain tumor through the catheter with a dimension of 0.9 mm, as shown in Figure 2c. [41,42] of interest with the dimension of 58 × 97 × 72 mm 3 is generated to fully enclose the brain tumor.
The final computational mesh consists of 82,944 uniform structured elements, obtained after the mesh-independence test [42]. The dimension of each element is 0.9 × 0.9 × 6.0 mm. Drugs are infused into the brain tumor through the catheter with a dimension of 0.9 mm, as shown in Figure 2c. Modelling results will be presented in the X-Z plane along the catheter track in the following sections, marked in green color.

Extraction of Tumor Heterogenous Properties
The quantitative analysis of DCE-MRI data is performed by computing the concentration of contrast agent from the signal intensity of the DCE-MR images on a voxel-byvoxel basis. This is based on the spoiled gradient recalled echo (SPGR)/FFE equation that is represented as [43] ( where (0) is the signal intensity when no contrast agent was given. (t) is the signal intensity of the contrast agent at a particular time point. Gd,ECS ( ) is the contrast agent concentration in the extracellular space (mmol • L −1 ). The longitudinal relaxivity ( 1 ) and transverse relaxivity ( 2 ) of Dotarem contrast agent in the plasma are 3.5 mmol −1 s −1 L and 4.9 mmol −1 s −1 L, respectively. The procedure developed in Ref [43] is adopted to calculate the pre-contrast 1 ( 10 ) using the three fast spin-echo (FSE) images, including T1weighted images ( ⁄ = 360 10 ⁄ ms), T2-weighted images ( ⁄ = 3500 90 ⁄ ms) and proton density (PD)-weighted images ( ⁄ = 3500 7.2 ⁄ ms). The rest imaging parameters for acquisition of T1-weighted, T2-weighted and PD-weighted images were the same as for T1-perfusion MR images, as summarised in Table 1.
The Leaky Tracer Kinetic Model (LTKM) [44] is applied to calculate the tissue perfusion parameters at each image voxel by fitting to the local contrast agent concentrations. Compared to the General Tracer Kinetic Model (GTKM), it is favored to analyze DCE-MRI data obtained in a short acquisition time [45,46], as the data used in this study (2 Modelling results will be presented in the X-Z plane along the catheter track in the following sections, marked in green color.

Extraction of Tumor Heterogenous Properties
The quantitative analysis of DCE-MRI data is performed by computing the concentration of contrast agent from the signal intensity of the DCE-MR images on a voxel-by-voxel basis. This is based on the spoiled gradient recalled echo (SPGR)/FFE equation that is represented as [43] I(t) where I(0) is the signal intensity when no contrast agent was given. I(t) is the signal intensity of the contrast agent at a particular time point. C Gd,ECS (t) is the contrast agent concentration in the extracellular space (mmol·L −1 ). The longitudinal relaxivity (R 1 ) and transverse relaxivity (R 2 ) of Dotarem contrast agent in the plasma are 3.5 mmol −1 s −1 L and 4.9 mmol −1 s −1 L, respectively. The procedure developed in Ref [43] is adopted to calculate the pre-contrast T 1 (T 10 ) using the three fast spin-echo (FSE) images, including T1-weighted images (TR/TE = 360/10 ms), T2-weighted images (TR/TE = 3500/90 ms) and proton density (PD)-weighted images (TR/TE = 3500/7.2 ms). The rest imaging parameters for acquisition of T1-weighted, T2-weighted and PD-weighted images were the same as for T1-perfusion MR images, as summarised in Table 1. The Leaky Tracer Kinetic Model (LTKM) [44] is applied to calculate the tissue perfusion parameters at each image voxel by fitting to the local contrast agent concentrations. Compared to the General Tracer Kinetic Model (GTKM), it is favored to analyze DCE-MRI data obtained in a short acquisition time [45,46], as the data used in this study (2 min). This is because GTKM requires relatively longer DCE-MRI data to stabilize the contrast agent's concentration and accurately determine the porosity. This is not always possible in the case of humans because of many reasons, such as prolonged scan time and distortion of the image caused by the patient's movement during scan time. The LTKM is expressed as where K trans is transvascular transport rate between plasma and extracellular space. λ tr is the volume transfer constant between plasma and leakage space. C Gd,BL as the arterial input function stands for the concentration of contrast agent in the blood, which can be obtained from the images following the procedure in Ref [47]. Readers are referred to the identified references for details of the image processing protocols.

Model Parameters
Given that the treatment period is much shorter than the tumor growth rate, the drug transport properties and tissue biological properties are treated as constants. Bevacizumab (BEV) and temozolomide (TMZ), as the typical antiangiogenic drug and cytotoxic drug, respectively, are selected in this study. The baseline values of model parameters are summarized in Tables 3 and 4 for tissues and drugs, respectively. Model parameters stranding for the heterogeneous tissue properties are extracted from DCE-MRI and mapped to each voxel of the 3D model geometry based on the coordinates. The MR images also consist of voxels present outside the brain, which belong neither to the tumor nor its surrounding tissue, i.e., exterior voxels corresponding to air. The model parameters are set as zero for these voxels [42,48].  [23] (S BL /V T ) b Baseline of blood vessel surface area per unit tissue volume m −1 2.0 × 10 4 [17] 7.0 × 10 3 [17] p BL Blood pressure Pa 4610 [49] 4610 [49] σ BL Osmotic reflection coefficient for proteins in blood − 0.82 [17] 0.91 [17] π BL Osmotic pressure of blood  The initial volume fraction of plasma ( BL,0 ) and extracellular space ( ECS,0 ) can be directly obtained from DCE-MR data as described in Section 2.2.2. Results of each image slice are given in Figure 3, where the tumor region presents a higher volume fraction of ECS than normal tissue. Since the ratio of CM / ICS is around 0.154 and 0.188 [22,65] for tumor cells and normal cells, respectively, the initial volume fraction of cell membrane ( CM,0 ) and intracellular space ( ICS,0 ) at each image voxel can be calculated by using Cancers 2022, 14, x FOR PEER REVIEW 9 of 33 Microvasculature surface becomes highly nonuniform as tumors grow. The ratio of trans (shown in Figure 3 for each slice) at each image voxel to the average value for the entire tumor region is then used to calculate the initial values for microvasculature density ( BL,0 T ⁄ ) by scaling with the standard value [68], ( BL T ⁄ ) b from the literature, as

Infusion Settings
Infusion settings including the infusion rate ( in ), infusion duration ( in ) and infusion concentration ( in ) are the factors that can be precisely controlled in clinical practice, directly determining the drug dose for administration. The infusion rate is usually kept in the range of 0.5~10.0 μL min ⁄ [69,70] to provide effective drug administration, while avoiding potential tissue damage. Moreover, it is found that the CED infusion catheter can be left indwelling for several days when the infusion rate is kept below 5.0 μL min ⁄ [7]. Since drugs were infused at the rate of 5.0 μL min ⁄ for 2~5 days in the clinical trials [69], the CED infusion is conducted at the rate of 3.0 μL min ⁄ in this study, lasting for 3 days. Given BEV and TMZ were infused at the concentrations of 7.  The antiangiogenic drugs, which are continuously infused into the brain tumor, can reduce the surface area of microvasculature over time by the scaling factor ϕ, as S BL /V T = ϕS BL,0 /V T . S BL is the surface area of the microvasculature wall, defined as S BL = 2πrl, where r and l are the mean radius and mean length of the blood vessels, respectively. Studies using microscopy imaging demonstrated that antiangiogenic drugs can not only successfully lead to the shrinkage of microvasculature, but also reduce the blood vessel length [66,67]. However, given there is no kinetic model to describe the relationship between the concentration of antiangiogenic drugs and the radius and length of blood vessels, it is further assumed that only the vessel radius reduces in response to the antiangiogenic drugs. Given the tissue volume remains constant, r = ϕr 0 . Moreover, since the volume fraction of plasma is defined as BL = V BL /V T = πr 2 l/V T , the antiangiogenic effect on the volume fraction of plasma ( BL ) can be expressed as The parameter, volume fraction of plasma ( BL ) is applied as the measure of microvascular density in this study. The value of ECS can then be updated using Equation (13) by treating the volume fractions of the cell membrane and intracellular space as constants.

Microvasculature Surface Area per Unit Tissue Volume (S BL /V T )
Microvasculature surface becomes highly nonuniform as tumors grow. The ratio of K trans (shown in Figure 3 for each slice) at each image voxel to the average value for the entire tumor region is then used to calculate the initial values for microvasculature density (S BL,0 /V T ) by scaling with the standard value [68], (S BL /V T ) b from the literature, as

Infusion Settings
Infusion settings including the infusion rate (R in ), infusion duration (T in ) and infusion concentration (C in ) are the factors that can be precisely controlled in clinical practice, directly determining the drug dose for administration. The infusion rate is usually kept in the range of 0.5 ∼ 10.0 µL/min [69,70] to provide effective drug administration, while avoiding potential tissue damage. Moreover, it is found that the CED infusion catheter can be left indwelling for several days when the infusion rate is kept below 5.0 µL/min [7]. Since drugs were infused at the rate of 5.0 µL/min for 2 ∼ 5 days in the clinical trials [69], the CED infusion is conducted at the rate of 3.0 µL/min in this study, lasting for 3 days. Given BEV and TMZ were infused at the concentrations of 7.26 × 10 −5 M [71] and 5.15 × 10 −3 M [72], respectively, in the preclinical experiments, the same settings are applied in the following simulations.

Numerical Methods
The governing equations are implemented into a finite volume method-based Computational Fluid Dynamics open-source code package, OpenFOAM, to generate numerical solutions. The transient term is discretized with the first-order Euler scheme. The Gauss linear, linear interpolation and Gauss linear upwind schemes are employed to discretize the gradient, Laplacian and divergence terms in the governing equations, respectively. Fluid pressure and velocity correction are linked by the pressure implicit splitting of operators (PISO) algorithm. The residual tolerances are set as 1.0 × 10 −6 for the interstitial fluid transport model and 1.0 × 10 −8 for the drug transport model, respectively, to control the modelling solution convergence. A fixed time step of 1.0 × 10 −3 s is selected after the time step-independence test. The governing equations of interstitial fluid flow are solved under the condition of no infusion in the first place to generate a steady-state solution, which represents the hydraulic environment in the tissues before the CED infusion takes place. This solution is then applied at time zero for simulating the interstitial fluid flow and drug transport upon CED infusion in a transient manner. The initial concentrations of BEV, liposomal TMZ and free TMZ are set as zero in the entire domain.
The numerical simulations are performed using a 64-bit Intel(R) Core i7-10700 processor (Clock speed: 2.90 GHz), eight cores with 32 GB RAM. The total computational time involved in solving the fluid flow and drug transport equations is approximately 1.5 h.

Boundary Conditions
The external surface of normal tissue is assumed to be fixed with zero gauge pressure and drug flux, as it is the source-sink-driven flow in tissue [41,42]. The continuous condition is employed at the tumor-normal tissue interface. The catheter wall is assumed to be rigid with zero slip or flux. The fixed infusion rate is specified at the infusion site, where the flux of free drugs is zero. The concentration of liposomes at the infusion site remains constant over the entire infusion duration for drug administration.

Qualification of Delivery Outcomes
The delivery outcomes of this combination therapy under different regimens are evaluated by the following qualification indexes from the perspectives of drug accumulation, drug distribution and treatment effectiveness.

Spatial-Averaged Concentration
The drug concentration directly reflects the amount of drugs accumulating in the tissue. It is determined by the aforementioned drug transport processes and varies throughout the brain tumor and its holding tissue. Spatial-averaged concentration (C FD,ECS,avg ) is applied to examine the drug accumulation in the entire tissue, which is expressed as where C FD ECS,i and V ECS,i are the local concentration of free drugs and the local tissue ECS volume, respectively. ECS,i is the volume fraction of local tissue ECS, and V T,i is local tissue volume. V ECS,T is the entire tissue ECS volume. The symbol 'i' refers to each control volume of the 3D computational mesh.

Location of Distance Course
The tumor microenvironment can vary considerably from the infusion site to the brain tumor periphery, presenting a heterogenous asymmetrical distribution. To evaluate this spatial change, the location of the distance course (ψ dis ) is applied, defined as in which ψ is the variable of interest. d i refers to the distance from the local tissue volume (V T,i ) to the infusion site, and dis is a given distance.

Distribution Non-Uniformity
The non-uniformity of drug distribution is represented by a dimensionless number, NUN, which is defined as NUN evaluates the spatial variation of the drug concentration. A lower value indicates a more uniform drug distribution.

Effective Distribution Volume
The drug effective distribution volume indexes the treatment effectiveness (V eff ) in which the local drug concentration is above the drug LD 90 , defined as

Model Validation
The modelling predicted drug delivery outcomes are compared with experiments in Figure 4. The model in Section 2.1 was applied, while treating brain tissue as a porous medium with homogenous properties; the corresponding values are summarized in Table 3. The drug diffusivity in brain tissue is set as 1.8 × 10 −11 m 2 /s [73], and its elimination rate is 1.0 × 10 −3 s −1 [73]. The modelling results agree with the experimental measurements [74], validating the model's accuracy in predicting the drug delivery outcomes.
The modelling predicted drug delivery outcomes are compared with experi Figure 4. The model in Section 2.1 was applied, while treating brain tissue as medium with homogenous properties; the corresponding values are summarized 3. The drug diffusivity in brain tissue is set as 1.8 × 10 −11 m 2 s ⁄ [73], and its eli rate is 1.0 × 10 −3 s −1 [73]. The modelling results agree with the experimental m ments [74], validating the model's accuracy in predicting the drug delivery outco

Baseline Study
The enhanced bulk flow of interstitial fluid is crucial for drug transport and lation in CED treatment. Due to antiangiogenesis, the dynamic change in tissue would reshape the hydraulic environment in the tumor and thereby alter the dr ery outcomes. The infusion catheter is positioned at the site where the microvas is densest. Results of simulations using Patient 1 data are represented below. Re Patient 2 show the same qualitative trends, which are included in Appendix A.

Antiangiogenic Effect
Antiangiogenic drugs are simultaneously infused with cytotoxic drugs into tumor in the combination therapy. Figure 5 represents the time course of BEV co tion in the entire tumor. BEV rapidly accumulates in the tumor on Day 1 upon th uous infusion. However, this increase slows down as time proceeds and eventua to a constant level on Day 3. This is due to the dynamic equilibrium between th term accounting for the drug supply by infusion and the sink term of drug elimi

Baseline Study
The enhanced bulk flow of interstitial fluid is crucial for drug transport and accumulation in CED treatment. Due to antiangiogenesis, the dynamic change in tissue structure would reshape the hydraulic environment in the tumor and thereby alter the drug delivery outcomes. The infusion catheter is positioned at the site where the microvasculature is densest. Results of simulations using Patient 1 data are represented below. Results for Patient 2 show the same qualitative trends, which are included in Appendix A.

Antiangiogenic Effect
Antiangiogenic drugs are simultaneously infused with cytotoxic drugs into the brain tumor in the combination therapy. Figure 5 represents the time course of BEV concentration in the entire tumor. BEV rapidly accumulates in the tumor on Day 1 upon the continuous infusion. However, this increase slows down as time proceeds and eventually tends to a constant level on Day 3. This is due to the dynamic equilibrium between the source term accounting for the drug supply by infusion and the sink term of drug elimination.
Cancers 2022, 14, x FOR PEER REVIEW 12 of Figure 5. The spatial-averaged concentration of BEV in the brain tumor as a function of time.
The spatial distribution of BEV on Day 3 is shown on a vertical plane along the cath eter track in Figure 6. The BEV concentration achieves its peak at the infusion site an decreases radially towards the domain boundary. Notably, although the infusion is highl directional, as pointed downward in Figure 6, the drugs can still transport along the cath eter track and accumulate posteriorly at the infusion site. This is related to the drug diffu  The spatial distribution of BEV on Day 3 is shown on a vertical plane along the catheter track in Figure 6. The BEV concentration achieves its peak at the infusion site and decreases radially towards the domain boundary. Notably, although the infusion is highly directional, as pointed downward in Figure 6, the drugs can still transport along the catheter track and accumulate posteriorly at the infusion site. This is related to the drug diffusive transport that is determined by the concentration gradient between the infusion site and the tissue. The spatial distribution of BEV on Day 3 is shown on a vertical plane along the catheter track in Figure 6. The BEV concentration achieves its peak at the infusion site and decreases radially towards the domain boundary. Notably, although the infusion is highly directional, as pointed downward in Figure 6, the drugs can still transport along the catheter track and accumulate posteriorly at the infusion site. This is related to the drug diffusive transport that is determined by the concentration gradient between the infusion site and the tissue. The tumor's biological properties before and after the combination therapy are compared on the same vertical plane in Figure 7. Results demonstrate the effectiveness of BEV in reducing the microvasculature surface area in unit tissue ( BL T ⁄ ) and microvascular density ( BL ). The antiangiogenic effect is not limited to the area in front of the catheter. Since BEV can transport deep in the tumor tissue as shown in Figure 6, microvasculature surface area and microvascular density around the catheter can also be reduced. This reduction can inhibit the drug loss by blood drainage, thereby retaining more drugs within the tumor for treatment. On the other hand, the volume fraction of extracellular space ( ECS ) can consequently be enlarged, which would decrease the tissue resistance to the transport of cytotoxic drugs into the deep tumor region. The tumor's biological properties before and after the combination therapy are compared on the same vertical plane in Figure 7. Results demonstrate the effectiveness of BEV in reducing the microvasculature surface area in unit tissue (S BL /V T ) and microvascular density ( BL ). The antiangiogenic effect is not limited to the area in front of the catheter. Since BEV can transport deep in the tumor tissue as shown in Figure 6, microvasculature surface area and microvascular density around the catheter can also be reduced. This reduction can inhibit the drug loss by blood drainage, thereby retaining more drugs within the tumor for treatment. On the other hand, the volume fraction of extracellular space ( ECS ) can consequently be enlarged, which would decrease the tissue resistance to the transport of cytotoxic drugs into the deep tumor region.
The antiangiogenic effects of BEV on tumor properties are further quantitatively evaluated in Figure 8. The upper panel shows that S BL /V T and microvascular density decreases exponentially over time. The relatively slow changes on Day 3 indicate that the further prolongation of BEV infusion results in minor modifications to the tissue microstructure. Furthermore, an opposite trend can be observed for the volume fraction of ECS. This is because the space the blood vessels release is occupied by the interstitial fluid, becoming the extracellular space.
The tumor properties are represented in the lower panel as a function of the distance from the infusion site. The values are calculated for discrete voxels along a sphere surface described by the same distance of d i . Comparisons show that BEV can successfully reduce the S BL /V T and microvascular density to a significantly low level in the entire tumor. Therefore, the volume fraction of ECS becomes higher throughout the brain tumor, easing the barrier for cytotoxic drugs to penetrate the tumor tissue. The antiangiogenic effects of BEV on tumor properties are further quantitatively evaluated in Figure 8. The upper panel shows that BL T ⁄ and microvascular density decreases exponentially over time. The relatively slow changes on Day 3 indicate that the further prolongation of BEV infusion results in minor modifications to the tissue microstructure. Furthermore, an opposite trend can be observed for the volume fraction of ECS. This is because the space the blood vessels release is occupied by the interstitial fluid, becoming the extracellular space.
The tumor properties are represented in the lower panel as a function of the distance from the infusion site. The values are calculated for discrete voxels along a sphere surface described by the same distance of . Comparisons show that BEV can successfully reduce the BL T ⁄ and microvascular density to a significantly low level in the entire tumor. Therefore, the volume fraction of ECS becomes higher throughout the brain tumor, easing the barrier for cytotoxic drugs to penetrate the tumor tissue. 3.2.2. Interstitial Fluid Flow Figure 9 compares the hydraulic environment before and after the combination therapy on the same vertical plane. Results show that IFP is higher in the tumor than in the surrounding normal tissue, consistent with the reported experiment finding [75]. Such high pressure is due to abnormal tumor properties, including high microvasculature sur-  Figure 9 compares the hydraulic environment before and after the combination therapy on the same vertical plane. Results show that IFP is higher in the tumor than in the surrounding normal tissue, consistent with the reported experiment finding [75]. Such high pressure is due to abnormal tumor properties, including high microvasculature surface area and blood vessel permeability, as shown in Table 3. However, as one moves towards the tumor periphery, IFP falls sharply at the interface between the tumor and normal tissue. More importantly, the CED infusion can override the original interstitial fluid flow to generate even higher pressure at the infusion site and raise the pressure in the entire tumor. Quantitative comparisons between the pre-and post-treatment hydraulic environments are given in Figure 10. IFP is uniform throughout the brain tumor before the combination treatment starts. This is because of the dynamic equilibrium reached between the blood pressure, osmotic pressure and IFP. In contrast, CED infusion would dramatically increase IFP around the infusion site to over 1.0 × 10 5 Pa. Although this pressure experiences a rapid drop in approximately 4 mm away from the infusion site, the pressure in the whole tumor can still be raised. Thus, IFP throughout the tumor gradually increases over time in response to the continuous CED infusion. Similar trends can be found in the distance courses of IFV. The velocity reaches its peak at the infusion site and quickly decreases with the distance; however, it remains relatively higher in the deep tumor tissue (as seen in the zoomed portion). The time course of intratumoral spatial-averaged IFV exhibits two-phase changes. It increases linearly in the first 4~5 h after the CED infusion starts, while the increase slows down as the infusion continues. Moreover, modelling results show that the fluid leakage decreases once the antiangiogenic drugs are infused. The combination therapy can effectively reduce the fluid leakage from the blood in the entire tumor, as indicated by its distance course shown in the lower panel. Relatively high IFV only occurs at the tumor-normal tissue interface due to the sharp fall of local IFP, which also promotes the outward convective transport of drugs from the tumor to its holding tissue. The velocity inside the tumor is low due to the even IFP distribution [17]. In contrast, CED infusion can accelerate the interstitial fluid flow in the entire brain tumor, with the most significant enhancement occurring at the infusion site. This rapid flow would improve convective drug transport for deep tissue penetration. Moreover, the comparison between the pre-and post-treatment denotes that the fluid leakage from blood can be significantly reduced in the combination therapy, due to the reduction of microvasculature density. This inhibited fluid leakage is beneficial for preventing the drug concentration from being diluted.

Interstitial Fluid Flow
Quantitative comparisons between the pre-and post-treatment hydraulic environments are given in Figure 10. IFP is uniform throughout the brain tumor before the combination treatment starts. This is because of the dynamic equilibrium reached between the blood pressure, osmotic pressure and IFP. In contrast, CED infusion would dramat-ically increase IFP around the infusion site to over 1.0 × 10 5 Pa. Although this pressure experiences a rapid drop in approximately 4 mm away from the infusion site, the pressure in the whole tumor can still be raised. Thus, IFP throughout the tumor gradually increases over time in response to the continuous CED infusion. Similar trends can be found in the distance courses of IFV. The velocity reaches its peak at the infusion site and quickly decreases with the distance; however, it remains relatively higher in the deep tumor tissue (as seen in the zoomed portion). The time course of intratumoral spatial-averaged IFV exhibits two-phase changes. It increases linearly in the first 4~5 h after the CED infusion starts, while the increase slows down as the infusion continues. Moreover, modelling results show that the fluid leakage decreases once the antiangiogenic drugs are infused. The combination therapy can effectively reduce the fluid leakage from the blood in the entire tumor, as indicated by its distance course shown in the lower panel.

Drug Concentration
The performance of the combination therapy using liposomal TMZ and BEV is evaluated by comparing to another three regimens under identical infusion settings. These include plain TMZ infusion, plain TMZ and BEV combined infusion, and liposomal TMZ infusion. Figure 11 compares the spatial distribution of free TMZ at different time points. The infusion of plain TMZ alone leads to minimal drug accumulation. This invisible drug concentration can be attributed to the fast drug elimination by blood drainage and bioreactions. Both the co-delivery with BEV and the application of liposomes can improve the delivery results. This is because the former inhibits blood drainage by reducing the microvasculature surface area, while the latter protects the drugs from unpreferred reactions. The most effective drug delivery can be obtained by co-delivering the liposomal TMZ with BEV. Moreover, the drug concentration presents similar distribution patterns since Day 2. This implies that the drug transport would reach a quasi-steady state after 48 h since the infusion starts; further prolongation of the infusion will contribute less to the distribution of drugs in the brain tumor.
Another set of simulations is performed based on the assumption of homogeneous tissue and drug properties. The microvascular density, tissue porosity and cell density are

Drug Concentration
The performance of the combination therapy using liposomal TMZ and BEV is evaluated by comparing to another three regimens under identical infusion settings. These include plain TMZ infusion, plain TMZ and BEV combined infusion, and liposomal TMZ infusion. Figure 11 compares the spatial distribution of free TMZ at different time points. The infusion of plain TMZ alone leads to minimal drug accumulation. This invisible drug concentration can be attributed to the fast drug elimination by blood drainage and bioreactions. Both the co-delivery with BEV and the application of liposomes can improve the delivery results. This is because the former inhibits blood drainage by reducing the microvasculature surface area, while the latter protects the drugs from unpreferred reactions. The most effective drug delivery can be obtained by co-delivering the liposomal TMZ with BEV. Moreover, the drug concentration presents similar distribution patterns since Day 2. This implies that the drug transport would reach a quasi-steady state after 48 h since the infusion starts; further prolongation of the infusion will contribute less to the distribution of drugs in the brain tumor.

comes.
The time courses of spatial-averaged free TMZ concentration are compared between different regimens in Figure 12a. The drug concentration increases sharply in the first 12 h and remains high over time. The maximum accumulation is achieved by co-delivery of liposomal TMZ with BEV, followed by the co-delivery of plain TMZ and BEV, liposomal TMZ infusion and plain TMZ infusion in sequence. A similar order can be found for the distance courses in Figure 12b. Furthermore, regardless of the dosing regimen, the drug concentrations reach their peak values at the infusion site and decrease rapidly with the increase in distance. The high TMZ concentrations can only be found in the area approximately 8 mm from the infusion site, indicating the delivery outcomes of CED are highly localized.  Another set of simulations is performed based on the assumption of homogeneous tissue and drug properties. The microvascular density, tissue porosity and cell density are set uniform throughout the tumor and its surrounding tissue with their spatiallyaveraged values, respectively. At the same time, the rest infusion settings are kept identical. Results are shown in Appendix B. Drugs present a more symmetrical profile as compared to Figure 11. Moreover, the drug concentrations are overpredicted under this idealized condition, particularly for the plain TMZ infusion. These findings further demonstrate the role of a heterogeneous intratumoral environment in determining drug delivery outcomes.
The time courses of spatial-averaged free TMZ concentration are compared between different regimens in Figure 12a. The drug concentration increases sharply in the first 12 h and remains high over time. The maximum accumulation is achieved by co-delivery of liposomal TMZ with BEV, followed by the co-delivery of plain TMZ and BEV, liposomal TMZ infusion and plain TMZ infusion in sequence. A similar order can be found for the distance courses in Figure 12b. Furthermore, regardless of the dosing regimen, the drug concentrations reach their peak values at the infusion site and decrease rapidly with the increase in distance. The high TMZ concentrations can only be found in the area approximately 8 mm from the infusion site, indicating the delivery outcomes of CED are highly localized.  The delivery outcomes of different regimens are further evaluated from the perspectives of drug distribution and treatment effectiveness in Figure 13. The results show that the direct infusion of plain TMZ leads to the most heterogeneous distribution in the brain tumor. Although the application of liposomes can homogenize drug distribution, this improvement is not as significant as co-delivery with BEV. The most uniform distribution can be found for the combined delivery of liposomal TMZ with BEV. The effective distribution volume varies distinctly between the regimens. Direct infusion of plain TMZ and liposomal TMZ fails to provide enough drug concentration to kill tumor cells effectively. In contrast, the treatment effectiveness can be improved by using antiangiogenic drugs. The liposomal TMZ and BEV combined delivery provides the most effective treatment. The delivery outcomes of different regimens are further evaluated from the perspectives of drug distribution and treatment effectiveness in Figure 13. The results show that the direct infusion of plain TMZ leads to the most heterogeneous distribution in the brain tumor. Although the application of liposomes can homogenize drug distribution, this improvement is not as significant as co-delivery with BEV. The most uniform distribution can be found for the combined delivery of liposomal TMZ with BEV. The effective distribution volume varies distinctly between the regimens. Direct infusion of plain TMZ and liposomal TMZ fails to provide enough drug concentration to kill tumor cells effectively. In contrast, the treatment effectiveness can be improved by using antiangiogenic drugs. The liposomal TMZ and BEV combined delivery provides the most effective treatment. Figure 13. Comparisons of (a) drug distribution non-uniformity, and (b) effective distribution volume between different regimens in the brain tumor on Day 3.

Effect of Heterogeneous Microvasculature
The distribution of microvasculature can vary considerably throughout a brain tumor, depending on the tumor location and growth stage. This heterogeneous tumor characteristic would directly influence the combination therapy where BEV is used to reduce the microvascular density. A statistical analysis shows that the microvascular density ( BL ) spans from 2.22 × 10 −14 to 9.66 × 10 −2 in this studied brain tumor. Therefore, three infusion sites with different BL are selected to examine the effects of microvasculature heterogeneity. The corresponding BL at these locations are 2.22 × 10 −14 (Case 1), 1.28 × 10 −2 (Case 2) and 9.66 × 10 −2 (Case 3), respectively. Figure 14 shows the time courses of free TMZ concentrations in the entire brain tumor Figure 13. Comparisons of (a) drug distribution non-uniformity, and (b) effective distribution volume between different regimens in the brain tumor on Day 3.

Effect of Heterogeneous Microvasculature
The distribution of microvasculature can vary considerably throughout a brain tumor, depending on the tumor location and growth stage. This heterogeneous tumor character-istic would directly influence the combination therapy where BEV is used to reduce the microvascular density. A statistical analysis shows that the microvascular density ( BL ) spans from 2.22 × 10 −14 to 9.66 × 10 −2 in this studied brain tumor. Therefore, three infusion sites with different BL are selected to examine the effects of microvasculature heterogeneity. The corresponding BL at these locations are 2.22 × 10 −14 (Case 1), 1.28 × 10 −2 (Case 2) and 9.66 × 10 −2 (Case 3), respectively. Figure 14 shows the time courses of free TMZ concentrations in the entire brain tumor for the infusion sites with different plasma volume fractions. The locations of the infusion site are shown in Figure 14a. After the treatment starts, TMZ begins to accumulate in the tumor ECS since drugs are continuously administrated. However, the accumulation rate varies. The concentration experiences a rapid increase in the first 24 h and remains at a relatively high level when the infusion site has the densest microvasculature. In contrast, infusing drugs into an area with fewer blood vessels results in a continuous increase of drug concentration over three days. Comparisons demonstrate that drug accumulation increases with the plasma volume fraction at the infusion site. The highest concentration can be achieved by placing the infusion catheter in the area with the highest BL .   Figure 15 compares the drug distribution and treatment effectiveness for the infusion sites with different plasma volume fractions. Results show that the distribution non-uniformity is inversely correlated to the local BL . Infusing drugs in the area with a higher BL enables drugs to transport into deeper tumor tissue, leading to a more uniform distribution. This is because the antiangiogenic drugs can effectively reduce the drug loss to the blood and enlarge the extracellular space for drugs to transport. Consequently, the treatment effectiveness increases with the local BL , since the deeper penetration allows the drugs to cover a larger tumor volume for effective cell killing.  Figure 15 compares the drug distribution and treatment effectiveness for the infusion sites with different plasma volume fractions. Results show that the distribution nonuniformity is inversely correlated to the local BL . Infusing drugs in the area with a higher BL enables drugs to transport into deeper tumor tissue, leading to a more uniform distribution. This is because the antiangiogenic drugs can effectively reduce the drug loss to the blood and enlarge the extracellular space for drugs to transport. Consequently, the treatment effectiveness increases with the local BL , since the deeper penetration allows the drugs to cover a larger tumor volume for effective cell killing. Figure 15 compares the drug distribution and treatment effectiveness for the infusion sites with different plasma volume fractions. Results show that the distribution non-uniformity is inversely correlated to the local BL . Infusing drugs in the area with a higher BL enables drugs to transport into deeper tumor tissue, leading to a more uniform distribution. This is because the antiangiogenic drugs can effectively reduce the drug loss to the blood and enlarge the extracellular space for drugs to transport. Consequently, the treatment effectiveness increases with the local BL , since the deeper penetration allows the drugs to cover a larger tumor volume for effective cell killing.

Effect of Heterogeneous Cell Density
Cell density can be highly heterogeneous in a tumor, particularly for the tumor at an advanced stage and with a large dimension. This tumor property, reflected by the volume fraction of ICS ( ICS ), affects the local drug cell uptake and the tumor ECS where the drugs transport. Given ICS is in the range of 5.808 × 10 −1 to 9.99 × 10 −1 in this examined tumor, three locations with different ICS are selected to examine its effects. The corresponding ICS at these three locations are 5.808 × 10 −1 (Case 1), 8.144 × 10 −1 (Case 2) and 9.99 × 10 −1 (Case 3), respectively.
The time courses of free TMZ concentration are compared for the infusion sites with different cell densities in Figure 16. Results show that the effectiveness of drug accumulation presents a negative relationship with the volume fraction of ICS at the infusion site. One possible reason is the low tissue porosity and microvascular density at the location where the cell density is high, resulting in less pronounced BEV-introduced improvement of drug transport. Cell density can be highly heterogeneous in a tumor, particularly for the tumor at an advanced stage and with a large dimension. This tumor property, reflected by the volume fraction of ICS (ϵ ICS ), affects the local drug cell uptake and the tumor ECS where the drugs transport. Given ϵ ICS is in the range of 5.808 × 10 −1 to 9.99 × 10 −1 in this examined tumor, three locations with different ϵ ICS are selected to examine its effects. The corresponding ϵ ICS at these three locations are 5.808 × 10 −1 (Case 1), 8.144 × 10 −1 (Case 2) and 9.99 × 10 −1 (Case 3), respectively.
The time courses of free TMZ concentration are compared for the infusion sites with different cell densities in Figure 16. Results show that the effectiveness of drug accumulation presents a negative relationship with the volume fraction of ICS at the infusion site. One possible reason is the low tissue porosity and microvascular density at the location where the cell density is high, resulting in less pronounced BEV-introduced improvement of drug transport. The impact of cell volume fraction on drug distribution and treatment effectiveness is represented in Figure 17. The distribution of free TMZ would be more non-uniform when the infusion catheter is placed at the area with a high ICS . Moreover, the effective The impact of cell volume fraction on drug distribution and treatment effectiveness is represented in Figure 17. The distribution of free TMZ would be more non-uniform when the infusion catheter is placed at the area with a high ICS . Moreover, the effective distribution volume decreases with ICS . These findings indicate that drug penetration is strongly limited, although infusing drugs to a location with higher cell density effectively kills the local tumor cells near the infusion site. This leads to a smaller tumor region with adequate drugs for treatment. Therefore, a better therapy could be achieved by infusing drugs to the location where the cell density is low.

Discussion
CED can effectively overcome the BBB by directly infusing drugs into the tumor. The rapid interstitial fluid flow improves the convective drug transport, enabling deep penetration into the tumor [5]. However, therapeutic effectiveness remains limited in clinical practice [6]. The intratumoral heterogeneous environment has been identified as one of the main limitations [7]. This modelling study demonstrates the effectiveness of co-delivery of cytotoxic drugs with antiangiogenic drugs for improving CED performance (Figures 11 and 12). On the one hand, the decreased microvasculature surface area and density reduce fluid leakage from the blood (Figure 9), preventing dilution of drug concentrations (Equation (7)). On the other hand, drug elimination by blood drainage can be inhibited (Equation (7)). Since the previous study showed the CED-infused drugs would transport into the blood and accumulate in other tissues and organs [76], this inhibition is beneficial to retain more cytotoxic drugs within the tumor for therapy. Moreover, the delivery outcomes can be further improved by using liposomal cytotoxic drugs ( Figure 13). This is because its low elimination rate allows the concentration of cytotoxic drugs to sustain at a relatively high level over time.
Given the heterogeneous characteristics of tumors, the location for positioning the catheter becomes critical. Treatment can be improved by infusing drugs at a site with denser microvasculature, as shown in Figure 15. As indicated in Figure 17, an inverse relationship is found between cell density and delivery outcomes; drug accumulation is more effective when the infusion catheter is placed at a site with a low cell density ( Figure  16).
CED infusion would result in a significant increase in pressure at the infusion site. It is predicted to be greater than 1.0 × 10 5 Pa, similar to the experimental measurements in the cat brains [77]. One must note that a brain tumor and its holding tissue can deform due to its soft nature in response to the CED infusion. Such high pressure can potentially cause tissue deformation and even tissue damage, consequently affecting drug delivery outcomes. Different models, including the poroelasticity model and hyper-viscoelasticity model, have been developed and applied to describe brain deformation [24,[78][79][80]. However, due to the lack of support to obtain heterogeneous mechanical properties of tumor tissue, the CED-induced tumor deformation and its influence on the delivery outcomes are not addressed in this study. This impact can be examined by performing a parameter analysis in the following study. To avoid large tissue deformation, the infusion rate must

Discussion
CED can effectively overcome the BBB by directly infusing drugs into the tumor. The rapid interstitial fluid flow improves the convective drug transport, enabling deep penetration into the tumor [5]. However, therapeutic effectiveness remains limited in clinical practice [6]. The intratumoral heterogeneous environment has been identified as one of the main limitations [7]. This modelling study demonstrates the effectiveness of co-delivery of cytotoxic drugs with antiangiogenic drugs for improving CED performance (Figures 11  and 12). On the one hand, the decreased microvasculature surface area and density reduce fluid leakage from the blood (Figure 9), preventing dilution of drug concentrations (Equation (7)). On the other hand, drug elimination by blood drainage can be inhibited (Equation (7)). Since the previous study showed the CED-infused drugs would transport into the blood and accumulate in other tissues and organs [76], this inhibition is beneficial to retain more cytotoxic drugs within the tumor for therapy. Moreover, the delivery outcomes can be further improved by using liposomal cytotoxic drugs ( Figure 13). This is because its low elimination rate allows the concentration of cytotoxic drugs to sustain at a relatively high level over time.
Given the heterogeneous characteristics of tumors, the location for positioning the catheter becomes critical. Treatment can be improved by infusing drugs at a site with denser microvasculature, as shown in Figure 15. As indicated in Figure 17, an inverse relationship is found between cell density and delivery outcomes; drug accumulation is more effective when the infusion catheter is placed at a site with a low cell density ( Figure 16).
CED infusion would result in a significant increase in pressure at the infusion site. It is predicted to be greater than 1.0 × 10 5 Pa, similar to the experimental measurements in the cat brains [77]. One must note that a brain tumor and its holding tissue can deform due to its soft nature in response to the CED infusion. Such high pressure can potentially cause tissue deformation and even tissue damage, consequently affecting drug delivery outcomes. Different models, including the poroelasticity model and hyper-viscoelasticity model, have been developed and applied to describe brain deformation [24,[78][79][80]. However, due to the lack of support to obtain heterogeneous mechanical properties of tumor tissue, the CED-induced tumor deformation and its influence on the delivery outcomes are not addressed in this study. This impact can be examined by performing a parameter analysis in the following study. To avoid large tissue deformation, the infusion rate must be well controlled. It is recommended to be kept below 10 µL/min to avoid tissue damage. The modelling results also denote that the increase in pressure is highly localized. The average pressure across the tumor increases gradually during the infusion period. A multiphysics model incorporating fluid mechanics, tissue mechanics and fluid-solid interaction will need to be developed for in-depth analysis. This combination therapy involves complex physiological and physicochemical processes that are determined by various tissue and drug properties. A comprehensive parameter analysis using a simplified idealized model geometry allows for evaluating the role of each factor to identify the most influential ones. Results can reveal the underlying mechanisms in CED and benefit the treatment design for better effectiveness.
It is worth noting that the delivery outcomes of CED treatment strongly depend on multiple factors, including the catheter shape, infusion direction and injection profile. Future studies can be carried out to examine their impact for optimization. Results can not only improve the design of CED medical devices, but also contribute to the development of treatment protocols and guidelines.
The model predictive power for simulating drug delivery has been validated in Figure 4 and in several past studies. The IFP in the solid tumor was predicted to be 1500 Pa [31], which was well within the experimental range of 1064~3990 Pa [81,82]. The modelling results showed that the IFV in brain parenchyma was 0.65 µm/s [23]; this was in agreement with experimental measurements [49,83]. For CED, the model-predicted penetration depths of nanoparticles were consistent with the measurements from the in vivo experiments [40]. Similar comparisons can also be found for small molecules when infused into the gel phantom [84] and animal brain [73]. However, one must note that the model applied in this study is developed to capture the key physiological and physicochemical processes involved in the CED treatment. The model parameters listed in Tables 3 and 4 are averaged and representative values from the literature. Therefore, the modelling results can only provide a qualitative trend of the delivery outcomes. Findings allow for identifying the importance of the examined influencing factors and determining the opportunities to improve the delivery outcomes. Medical imaging provides a non-invasive solution to obtain a realistic, patient-specific intratumoral environment, such as interstitial fluid flow and transport [85] and tissue anisotropy [35]. The image-derived information can be applied to improve the modelling accuracy.
This study involves several assumptions and limitations.
(a) The insertion of a catheter into the brain tissue has the potential to result in trauma and oedema. The fluid generated from damaged cells and the enhanced leakage from damaged blood vessels [86] can alter the interstitial fluid flow and cause swelling of brain tissue. However, there is a lack of accurate models to describe such a physiochemical process. Moreover, since oedema eventually disappears and the tissue properties can get back to normal levels [87,88], the effect of oedema is not considered. The mathematical model needs to be developed with support from in vivo experiments to describe this process. (b) The antiangiogenic effect is reflected by the reduction of blood vessel diameter. Although blood vessels shirk significantly in response to the antiangiogenic drugs [66], both the vessel diameter and length can change [67]. Disappointingly, there is no such model to define the relationship between the changes in morphological characteristics of microvasculature and the antiangiogenic drug concentration. For improvement, results from microscale studies on the dynamic changes of blood vessels in antiangiogenic treatment will be needed for model development. In addition, modelling on the capillary level can be applied to predict the dynamic response of the microvas-cular network to the antiangiogenic drugs [89,90], which would shed light on the development of the combination therapy. (c) Drug diffusivity can also be location-dependent in a brain tumor, subject to the local in vivo environment. This variation of drug property can also affect drug transport and accumulation. However, the diffusivity of each drug is assumed to be uniform in the brain tumor in this study due to the lack of relevant information that can be extracted from the applied medical images. This assumption can be relaxed by using diffusion-weighted MR images [91], where the signal intensity is related to the diffusivity of water molecules. (d) Liposome cell uptake is ignored in this study owing to its complex nature. This process is controlled by several factors including particle size [92], ligands [93] and/or energy consumption [94,95]. A bespoken model needs to be developed to describe this process when focusing on a specific type of liposome. (e) Only one catheter is placed in this study to examine the effectiveness of combination therapy in treating heterogeneous brain tumors. However, it is important to note that the cytotoxic drugs can penetrate a short distance from the infusion site; it is approximately 8 mm, as shown in Figure 12. This limited penetration depth results in highly localized cell killing, which is less effective for treating large tumors. In practice, it becomes more common to use multiple catheters [96] or a catheter with multiple injection-ports [97] simultaneously to enlarge the coverage of the drug to the tumor for better treatment outcomes. In this regard, the arrangement of catheters, the infusion directions and locations will be the key and worthy of in-depth study in the future. (f) The infusate administrated into the brain tumor is able to push the soft tissue back to form backflow. Drugs in the backflow would transport along the track of the infusion catheter to the normal tissue, reducing the treatment effectiveness. The formation of backflow depends on the infusion rate, tumor location and more importantly, the local tissue properties. Given the difficulty to obtain the heterogeneous tissue mechanical properties from the applied medical images, this phenomenon is not simulated in this study. This limitation can be overcome by incorporating the tissue mechanics model, solid-fluid interaction and corresponding tissue properties [27,98]. (g) In the current study, the anti-angiogenetic and cytotoxic drugs were infused simultaneously inside the brain tumor due to the lack of clinical information regarding their infusion timings. Future studies in this direction should attempt to optimize the infusion schedule between antiangiogenetic and cytotoxic treatment to maximize the therapeutic output.

Conclusions
The combined CED of liposomal cytotoxic drugs and antiangiogenic drugs to a heterogeneous brain tumor has been studied by image-based mathematical modelling. Results show that the co-infusion of antiangiogenic drugs can effectively improve the delivery outcomes of cytotoxic drugs. The reduced microvascular density inhibits the fluid exchange between the blood and tissue as well as the drug loss due to blood drainage. Moreover, the location of the infusion site plays an important role in determining the performance of this combination therapy. Better treatment can be achieved by infusing drugs at denser microvasculature sites. Infusing the drug in a place with high cell density risks reducing its effectiveness against the entire tumor. Results obtained from this study can deepen the understanding of chemotherapy using CED. Please note that this study is based on the DCE-MRI data of two patients, thus further studies using a larger patient cohort will be needed to validate the results and to provide a reference for improving the treatment design. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, upon reasonable request. posomal TMZ infusion, and liposomal TMZ and BEV combination infusion. It can be found that the infusion of free TMZ yields the worst drug profile, followed by liposomal TMZ infusion. The co-delivery with BEV can improve drug distribution; the best delivery outcomes are achieved by infusing BEV and liposomal TMZ. The time course and distance course of spatially averaged concentration of free TMZ in the brain tumor are compared between the abovementioned four delivery scenarios, in Figure A4. Results demonstrate the advantages of combination therapy using BEV and liposomal TMZ in improving drug accumulation and penetration.
The distribution heterogeneity of free TMZ is compared between the four delivery scenarios in Figure A5a. The combination delivery of BEV and liposomal TMZ leads to the most homogeneous distribution pattern; whereas drugs would concentrate in a limited region when free TMZ is mono-infused. Similarly, to the results in Figure 13, the largest effective distribution volume can be found for the co-infusion of BEV and liposomal TMZ, as denoted in Figure A5b. The only difference is the infusion of liposomal TMZ that can also provide effective cell killing, which is not shown in Figure 13. However, this effective distribution volume is still much lower than that of combination therapy. The time course and distance course of spatially averaged concentration of free TMZ in the brain tumor are compared between the abovementioned four delivery scenarios, in Figure A4. Results demonstrate the advantages of combination therapy using BEV and liposomal TMZ in improving drug accumulation and penetration.
The distribution heterogeneity of free TMZ is compared between the four delivery scenarios in Figure A5a. The combination delivery of BEV and liposomal TMZ leads to the most homogeneous distribution pattern; whereas drugs would concentrate in a limited region when free TMZ is mono-infused. Similarly, to the results in Figure 13, the largest effective distribution volume can be found for the co-infusion of BEV and liposomal TMZ, as denoted in Figure A5b. The only difference is the infusion of liposomal TMZ that can also provide effective cell killing, which is not shown in Figure 13. However, this effective distribution volume is still much lower than that of combination therapy.      Figure A5. Comparisons of (a) drug distribution non-uniformity, and (b) effective distribution v ume between different regimens in the brain tumor. The results are taken on Day 3.

Appendix A.3. Effect of Microvascular Density
Three infusion sites with different BL are selected in the second brain tumor, 4.6 × 10 −4 (Case 1), 0.01 (Case 2) and 0.04 (Case 3), respectively. The locations of the infusion sites are illustrated in Figure A6a. Results in Figure A6b denote that the accum lation of free TMZ is positively related to the local microvascular density, which is co sistent with the findings obtained from Figure 14. The delivery outcomes are evaluated in terms of NUN and eff in Figure A7. Sim lations show the free TMZ would present a more homogenous distribution pattern cover a larger tumor area for effective therapy when the infusion catheter is located at higher microvascular density site. This agrees with the results in Figure 15. Figure A5. Comparisons of (a) drug distribution non-uniformity, and (b) effective distribution volume between different regimens in the brain tumor. The results are taken on Day 3.

Appendix A.3. Effect of Microvascular Density
Three infusion sites with different BL are selected in the second brain tumor, as 4.6 × 10 −4 (Case 1), 0.01 (Case 2) and 0.04 (Case 3), respectively. The locations of these infusion sites are illustrated in Figure A6a. Results in Figure A6b denote that the accumulation of free TMZ is positively related to the local microvascular density, which is consistent with the findings obtained from Figure 14.

Appendix A.3. Effect of Microvascular Density
Three infusion sites with different BL are selected in the second brain tumor, as 4.6 × 10 −4 (Case 1), 0.01 (Case 2) and 0.04 (Case 3), respectively. The locations of these infusion sites are illustrated in Figure A6a. Results in Figure A6b denote that the accumulation of free TMZ is positively related to the local microvascular density, which is consistent with the findings obtained from Figure 14. The delivery outcomes are evaluated in terms of NUN and eff in Figure A7. Simulations show the free TMZ would present a more homogenous distribution pattern to cover a larger tumor area for effective therapy when the infusion catheter is located at a higher microvascular density site. This agrees with the results in Figure 15. The delivery outcomes are evaluated in terms of NUN and V eff in Figure A7. Simulations show the free TMZ would present a more homogenous distribution pattern to cover a larger tumor area for effective therapy when the infusion catheter is located at a higher microvascular density site. This agrees with the results in Figure 15.

Appendix A.4. Effect of Cell Density
Three infusion sites with different ICS are selected in the second brain tumor, 0.878 (Case 1), 0.9603 (Case 2) and 0.9872 (Case 3), respectively. The places where t infusion catheter is located in the brain tumor are shown in Figure A8a. Results in Figu A8b denote that the accumulation of free TMZ is inversely correlated to the local cell de sity at the infusion site. The evaluation of drug distribution and potential treatment effectiveness of the thr cases is shown in Figure A9. It can be found that infusing drugs at the location with t highest cell density results in the most non-uniform distribution and the worst cell killin These are consistent with the findings shown in Figure 17.

Appendix A.4. Effect of Cell Density
Three infusion sites with different ICS are selected in the second brain tumor, as 0.878 (Case 1), 0.9603 (Case 2) and 0.9872 (Case 3), respectively. The places where the infusion catheter is located in the brain tumor are shown in Figure A8a. Results in Figure A8b denote that the accumulation of free TMZ is inversely correlated to the local cell density at the infusion site.

Appendix A.4. Effect of Cell Density
Three infusion sites with different ICS are selected in the second brain tumor, as 0.878 (Case 1), 0.9603 (Case 2) and 0.9872 (Case 3), respectively. The places where the infusion catheter is located in the brain tumor are shown in Figure A8a. Results in Figure  A8b denote that the accumulation of free TMZ is inversely correlated to the local cell density at the infusion site. The evaluation of drug distribution and potential treatment effectiveness of the three cases is shown in Figure A9. It can be found that infusing drugs at the location with the highest cell density results in the most non-uniform distribution and the worst cell killing. These are consistent with the findings shown in Figure 17. The evaluation of drug distribution and potential treatment effectiveness of the three cases is shown in Figure A9. It can be found that infusing drugs at the location with the highest cell density results in the most non-uniform distribution and the worst cell killing. These are consistent with the findings shown in Figure 17.

Appendix B. Impact of the Assumption of Homogeneous Tissue Properties
In this study, it is assumed that the tissue properties are homogeneously distribut in the brain tumor and normal tissue, respectively. These properties include the microva cular density, tissue porosity, cell density and microvasculature surface area per tiss volume. The modelling is based on Patient 1′s data, and the spatially averaged values the abovementioned properties are applied. The rest of the properties' values are the sam as listed in Tables 3 and 4. Figure A10 represents the free TMZ distribution under the four delivery scenarios different time points. The drugs present a more symmetrical distribution compared to t results in Figure 10. Furthermore, the drug concentration is also overestimated. This highly associated with the non-uniformly distributed microvasculature, which can ra idly eliminate the drugs before they penetrate deep into the tissue. This compariso demonstrates the impact of the tumor heterogeneous environment on drug transport.

Appendix B. Impact of the Assumption of Homogeneous Tissue Properties
In this study, it is assumed that the tissue properties are homogeneously distributed in the brain tumor and normal tissue, respectively. These properties include the microvascular density, tissue porosity, cell density and microvasculature surface area per tissue volume. The modelling is based on Patient 1 s data, and the spatially averaged values of the abovementioned properties are applied. The rest of the properties' values are the same as listed in Tables 3 and 4. Figure A10 represents the free TMZ distribution under the four delivery scenarios at different time points. The drugs present a more symmetrical distribution compared to the results in Figure 10. Furthermore, the drug concentration is also overestimated. This is highly associated with the non-uniformly distributed microvasculature, which can rapidly eliminate the drugs before they penetrate deep into the tissue. This comparison demonstrates the impact of the tumor heterogeneous environment on drug transport.

Appendix B. Impact of the Assumption of Homogeneous Tissue Properties
In this study, it is assumed that the tissue properties are homogeneously distributed in the brain tumor and normal tissue, respectively. These properties include the microvas cular density, tissue porosity, cell density and microvasculature surface area per tissu volume. The modelling is based on Patient 1′s data, and the spatially averaged values o the abovementioned properties are applied. The rest of the properties' values are the sam as listed in Tables 3 and 4. Figure A10 represents the free TMZ distribution under the four delivery scenarios a different time points. The drugs present a more symmetrical distribution compared to th results in Figure 10. Furthermore, the drug concentration is also overestimated. This i highly associated with the non-uniformly distributed microvasculature, which can rap idly eliminate the drugs before they penetrate deep into the tissue. This comparison demonstrates the impact of the tumor heterogeneous environment on drug transport. Figure A10. Spatial distribution of released free temozolomide at different time points. Figure A10. Spatial distribution of released free temozolomide at different time points.