Modelling of combination therapy using implantable anticancer drug delivery with thermal ablation in solid tumor

Local implantable drug delivery system (IDDS) can be used as an effective adjunctive therapy for solid tumor following thermal ablation for destroying the residual cancer cells and preventing the tumor recurrence. In this paper, we develop comprehensive mathematical pharmacokinetic/pharmacodynamic (PK/PD) models for combination therapy using implantable drug delivery system following thermal ablation inside solid tumors with the help of molecular communication paradigm. In this model, doxorubicin (DOX)-loaded implant (act as a transmitter) is assumed to be inserted inside solid tumor (acts as a channel) after thermal ablation. Using this model, we can predict the extracellular and intracellular concentration of both free and bound drugs. Also, Impact of the anticancer drug on both cancer and normal cells is evaluated using a pharmacodynamic (PD) model that depends on both the spatiotemporal intracellular concentration as well as characteristics of anticancer drug and cells. Accuracy and validity of the proposed drug transport model is verified with published experimental data in the literature. The results show that this combination therapy results in high therapeutic efficacy with negligible toxicity effect on the normal tissue. The proposed model can help in optimize development of this combination treatment for solid tumors, particularly, the design parameters of the implant.


Modelling of combination therapy using implantable anticancer drug delivery with thermal ablation in solid tumor
Muneer Al-Zu'bi & Ananda Mohan * Local implantable drug delivery system (IDDS) can be used as an effective adjunctive therapy for solid tumor following thermal ablation for destroying the residual cancer cells and preventing the tumor recurrence. In this paper, we develop comprehensive mathematical pharmacokinetic/ pharmacodynamic (PK/PD) models for combination therapy using implantable drug delivery system following thermal ablation inside solid tumors with the help of molecular communication paradigm. Cancer is one of the most dangerous and deadliest diseases that cause deaths of millions of people around the world each year. More than 85% of human cancers appear as solid tumors 1 . Minimally invasive and imageguided thermal ablation techniques such as radiofrequency ablation (RFA) and high-intensity focused ultrasound (HIFU) ablation have been used for local treatment of malignant solid tumors as an alternative for systemic chemotherapy and surgical resection 2 . In thermal ablation, the tissue temperature rises above 50 °C, which is enough to destroy tumor tissue by inducing coagulation necrosis 2 . Thus, thermal ablation will destroy the cellular and vasculature structures, which results in a new ablated tissue with different characteristics. However, the main limitation and challenge of thermal ablation is the risk of local recurrence and residual tumor after the thermal ablation, particularly at the tumor periphery [3][4][5] . The experimental studies showed that the combination therapy through local release of anticancer drugs from a miniaturized implant in solid tumors following thermal ablation can result in a better therapeutic efficacy by destroying the residual cancer cells and preventing tumor recurrence [6][7][8] .

In this model, doxorubicin (DOX)-loaded implant (act as a transmitter) is assumed to be inserted inside solid tumor (acts as a channel) after thermal ablation. Using this model, we can predict the extracellular and intracellular concentration of both free and bound drugs. Also, Impact of the anticancer drug on both cancer and normal cells is evaluated using a pharmacodynamic (PD) model that depends on both the spatiotemporal intracellular concentration as well as characteristics of anticancer drug and cells. Accuracy and validity of the proposed drug transport model is verified with
Anticancer drug distribution and fluid flow within solid tumors are essential factors that affect the clinical efficiency of anticancer therapies. In addition to that, pharmacokinetic processes, including drug efflux/influx into cells, drug binding/unbinding with interstitial proteins, perfusion into blood capillaries, and biodegradation, have a significant impact on the therapeutic outcomes. Moreover, one of the most important challenges in the development of new drug delivery systems (DDSs) is ensuring that the optimal amount of drug is achieved in tumor versus normal tissues to avoid toxicity in the healthy tissues. Dual-release implants have been clinically approved and currently used for cancer treatments in the clinical trials and research, e.g., Gliadel polymer implants 9 . The implant releases the anticancer drug over two phases, namely, burst and sustained releases. For example, in-situ forming implants (ISFIs) have a dual-release pattern with a large undesirable burst release, which may cause major toxicity problems and consume the loaded drug in the implants rapidly 10 . Designing this type of implants with a minimum initial burst release becomes an attractive challenge, and one of the key issues in the design of ISFIs 11 . Moreover, some drug-loaded implants can be designed in a controlled way to www.nature.com/scientificreports/ medium surrounded by normal tissue, which acts as a molecular communication channel. At the target sites, the intracellular concentration of the anticancer drug (DOX) should reach a minimum threshold to kill the cancer cells. In the MC paradigm, this can be considered as a reception mechanism where the intracellular concentration is the received signal, while the death of cancer cells is the output response.
In this model, a millimeter-scale dual-release implant, loaded with anticancer DOX drug, is assumed to be inserted inside a solid tumor to releases DOX anticancer agents. Here, we consider two solid tumor models, namely, thermally ablated and non-ablated tumors, surrounded by normal tissues. We consider the impact of the following factors on the drug transport process in tumor and surrounding normal tissue: interstitial fluid pressure and velocity, binding of DOX with albumin-proteins in the interstitial extracellular space, cellular influx/efflux of DOX across the cellular membranes, elimination of DOX into the blood and lymphatic microvessels. Impacts of all the above-mentioned pharmacokinetic processes are included in the drug transport model for predicting the extracellular and intracellular concentrations of both free and bound DOX. Furthermore, this model enables estimate the toxicity of DOX on tumor cells and surrounding healthy tissue. The impact of DOX on the cancer cells is evaluated using a pharmacodynamic model that depends on the spatiotemporal intracellular concentration of DOX as well as on the characteristics of both the DOX and tumor cells. Moreover, the concentration of DOX in normal tissue is evaluated, which can be used for toxicity assessment. Accuracy and validity of our proposed model are verified and compared with the published experimental data in the literature, assuming the impact of the various pharmacokinetic parameters are combined in the apparent diffusivity and apparent elimination constant. To the best of our knowledge, this work is the first comprehensive model available in the literature that simultaneously captures and addresses the anticancer drug transport, pharmacokinetics, and pharmacodynamics using local dual-release drug implants in malignant solid tumors following thermal ablation.

Results and discussion
In this study, the governing mathematical equations which describe the proposed model are discretized in space with the finite element method using the commercial software package COMSOL Multiphysics 5.3. In COM-SOL Multiphysics, we solve the interstitial fluid flow and drug transport models together with the tumor cell density model. The steady-state solutions of the interstitial velocity and pressure fields are applied to the drug transport model. Both the fluid flow and the drug transport models are solved under relative tolerance of 10 -6 and absolute tolerance 10 -7 . The numerical computation is run to examine the model over a timeframe of 4 days (96 h), assuming the initial time is the time of insertion of the implant in the tumor. In this model, the tumor characteristics and the drug parameters are taken from the published experimental works and other studies in the literature. The parameters used in this study are given in Tables 2, 3, 4   www.nature.com/scientificreports/ As shown in Fig. 3, the average and maximum interstitial fluid pressure (IFP) in the tumor are equal to 1466 Pa and 1533 Pa, respectively. However, the average pressure in the normal tissue is equal to 50 Pa, which is significantly lower than the pressure in the tumor. Furthermore, there is a sharp pressure gradient (and consequently, high-velocity field) at the interface boundary between the tumor and normal tissues, as shown in Fig. 3. This result is not surprising due to the lake of the lymph vasculature in the tumor compared to the normal tissue. The trend of the results agrees well with the previous studies in the literature 1,27,28 .
A group of researchers conducted experiments for measuring the DOX concentration with the distance from polymer implants placed inside thermally ablated liver tumor 8,20,21,29 . In Fig. 4, we verify the accuracy and validity of our drug transport model by comparing the results with the published experimental data 21 . The impact of the various processes in the tumor, including the binding effect, is given in terms of the apparent elimination rate constant and apparent diffusivity 20 . The parameters used in this comparison are chosen to be the same as that used in the experiment 21 . The apparent diffusivity of DOX in the liver tumor is given as D = 50 µm 2 /s while it varies within the thermally ablated tumor; thus, we use an average value 20 of 78.2 µm 2 /s. The apparent elimination rate constant of DOX in the liver tumor is γ = 0.58 × 10 -4 s −1 , and it is negligible in the ablated tumor within the first 4 days 20 , i.e., γ = 0 s −1 . As expected, the measured concentration shows a decreasing trend with the radial distance from the implant. The results obtained using our numerical COMSOL model agree well with the results extracted from the published experimental data. This indicates the accuracy and validity of the drug release and transport model, which represents the main part of the proposed model in this paper.
In this study, the concentration distribution profiles for both 80% and 90% ablation have a similar trend with the time and distance but with slightly different amplitudes. Therefore, we do not show all the results in this paper to eliminate the redundancy. Figure 5 shows the spatial-mean extracellular concentration profiles of free-DOX and bound-DOX in the risk region of 90% ablated tumor. The implant with a higher release rate leads  www.nature.com/scientificreports/ to a larger peak concentration and lower peak time. As the release rate increases, the amount of the released DOX will increase rapidly, and thus the peak concentration will reach a higher value within a shorter time. The increasing and decreasing rates of the extracellular concentration become sharper as the release rate constant increases. The decay in concentration after it reaches the peak value is due to a reduction in the released drug from the implant and elimination through blood vessels and cellular uptake. The bound-DOX concentration has a similar trend as the free-DOX concentration for various sustained release rate constants but with approximately three-fold higher amplitude. Figure 6 shows the extracellular free-DOX concentration versus the distance from the implant/tissue interface in ablated and non-ablated solid tumors. The extracellular DOX concentration in a solid tumor without applying thermal ablation has a smaller amplitude than the concentration in an ablated tumor. This happens because in the non-ablated tumor, the cells and blood vasculature structures, which cover a large volume of the solid tumor, have a high impact on the elimination of the DOX through the cellular uptake and the blood microvessels. www.nature.com/scientificreports/ Furthermore, the drug can penetrate a larger distance and cover a larger volume in the ablated tumor compared to the non-ablated tumor. Also, tumor with larger ablation radius (e.g., 90%) shows higher DOX concentration and larger penetration compared to smaller ablation radius (e.g., 80%) on the fourth day as shown in Fig. 6. The observed impact of the thermal ablation on the drug distribution in tumors agrees with the experimental data in the literature 6,8,21 .
The intracellular concentration of free-DOX in the risk region of 90% ablated tumor is shown in Fig. 7. The DOX intracellular concentration follows a similar trend as the extracellular concentration because it highly depends on the extracellular DOX levels. The higher peak amplitude of the intracellular concentration can be achieved using a faster release implant. As shown in Fig. 7b, the intracellular concentration decreases as the distance increases from the inner boundary of the risk region. This happens because the extracellular DOX concentration decreases with the distance, and it has a direct influence on the intracellular uptake and, consequently, on the intracellular DOX concentration. Figure 8 shows the tumor cell density in the risk region with the radial distance from the inner boundary of the risk region of 80% ablated tumors at various times following the insertion of the implant. Tumor cell density shows heterogeneous distribution along the radial direction, where it increases with the distance with minimum density appears near the inner boundary of the risk region. This can also be seen in the color map of the spatial distribution of tumor cell density in Fig. 10. Moreover, there is a significant reduction in tumor cell density over time following the insertion of the implant. Furthermore, we found that the implant with the release rate constants k s = {5, 25} × 10 -6 s −1 will almost have the same therapeutic effect on the last day of treatment. However, using a smaller release rate will consume less amount of the drug with minimum toxicity on the normal tissue. Therefore, design the implant with an optimal release rate is very important to get a high therapeutic efficacy.
The final therapeutic outcomes of the combination therapy using the DOX-loaded implant and RFA can be obtained from the tumor survival curves, as shown in Fig. 9. We can see that using the implant alone for tumor treatment without RFA leads to a negligible impact on the tumor cell density, even with a high release rate, i.e.,   Fig. 10. Thus, the combination treatment using the implant following RFA will result in high therapeutic efficacy in destroying the residual tumor cells compared to a single therapy approach. Moreover, the release rate constants, k s = {5, 25} × 10 −6 s −1 , show a similar therapeutic effect on the last day (at t = 96 h) with lower tumor survival compared to an implant with k s = 1 × 10 −6 s −1 .
In general, the IDDSs have negligible toxicity on the surrounding healthy tissue and do not cause systemic toxicity. However, the amount of DOX which may reach to the normal tissue should be minimized to reduce the toxicity risks. The DOX level in the normal tissue can be used as a metric to predict the toxicity. In this study, the peak DOX concentration in normal tissue under all the examined release rates, as listed in Table 1, is found to be lower than the half-maximal inhibitory concentration 27 (IC 50 = 4.13 × 10 −5 kg/m 3 ). The peak extracellular concentration of free and bound-DOX in the solid tumor without RFA has very small and negligible values. In the non-ablated solid tumor, the released DOX from the implant will be affected by the cellular uptake and the

Conclusions
In this paper, we propose comprehensive mathematical and computational models for the transport of anticancer drug following the insertion of a dual-release implant in thermally ablated solid tumor with the help of the molecular communication paradigm. We predict the extracellular and intracellular concentrations of both free and bound-DOX in the various regions of thermally ablated solid tumor. This model includes the impact of the various pharmacokinetic processes such as binding of DOX to proteins, cellular influx/efflux, and elimination into the vasculature system. Also, we investigate the impact of the pressure and velocity of interstitial fluid on DOX transport in tumor and surrounding healthy tissue. Accuracy and validity of the proposed transport model is verified with the published experimental data assuming that the various pharmacokinetic processes are combined in the apparent diffusivity and elimination constant. Moreover, we examine the impact of the anticancer drug on tumor cell density, which shows a significant reduction in cell density over time. The combination therapy using the implantable drug delivery following thermal ablation results in high therapeutic efficacy. We found that the anticancer drug does not lead to toxicity effect on the normal tissue. The proposed model can help to optimize the development of the combination technique for treating solid tumors. Thus, we can reduce the clinical trials and the number of animals in biomedical research to save time and reduce cost. One of the limitations in our model is ignoring the impact of the implant biodegradation and drug interactions with the tissue surrounding the implant on the drug release process. Another limitation is using average diffusivity within the ablated tumor. However, the diffusivity varies with the radial distance in the ablated tumor. We will improve our model to overcome these limitations in future works.

Method
Mathematical model. Solid tumors are heterogeneous environments due to spatial heterogeneity of the tumor vasculature and the cells. However, due to the unavailability of experimental heterogeneity data of solid tumors and to simplify the analysis, the solid tumors are widely treated in the literature as spatially homogeneous media 27,28,[30][31][32] . Thus, we do not discriminate between the necrotic and viable tumor regions. Moreover, we assume that the growth timescale of the tumor and normal tissues is much longer than the timescale of the transport phenomena and the observation time window. Thus, it would be reasonable to assume that the system's physiological parameters to be time-independent 27 . Incomplete radiofrequency ablation (RFA) will create an ablated zone, with no viable cancer cells, surrounded by a tumor rim (risk region) that shows a high density of the viable malignant cells, as shown in Fig. 11.
Tumor microenvironment modelling. Interstitial fluid transport. The tumor and surrounding tissue can be treated as porous media because the length scale of the intercapillary distances is much smaller than the tumor radius 27,30,33,34 . Thus, the variations over the microscopic length scales can be averaged out, and the interstitial fluid flow (IFF) is defined by coupling mass and momentum conservation equations. For incompressible Newtonian fluid flow through a porous medium, the momentum conservation equation (Navier-Stokes equation) is simplified to Darcy's law at a steady-state, which is quite applicable to the analysis of the interstitial fluid flow 33 . Darcy's law is used to account for the convective contribution of the interstitial fluid through porous media. Darcy law is derived for incompressible Newtonian fluid with neglecting the following factors: the divergence of the velocity, the inertial force, and the friction within the fluid and between the fluid and solid phases. Similar to other works in the literature 35,36 , the fluid flow in the tumor tissue with thermal ablation can be modelled using Darcy law.
where v i is the interstitial velocity field (IVF) in (m/s), P i is the interstitial fluid pressure (IFP) in (Pa), k i is the hydraulic conductivity of the interstitial fluid in (m 2 /(Pa s)).
The mass continuity (balance) equation for an incompressible fluid in the porous media with source and sink of mass is given as follows 27 . www.nature.com/scientificreports/ where ρ i is the interstitial fluid density in (kg/m 3 ), φ vi is the mass fluid source term which represents the fluid flow rate per unit volume of tissue from the blood vessels into the interstitial space in (1/s), and φ li is the lymphatic drainage (sink) term that represents the fluid flow rate per unit volume of tissue from the interstitial space into the lymph vessels in (1/s). The Eq. (2) is applicable to both normal and cancerous biological tissues. The mass fluid source term φ vi is governed by Starling's law 27,30 as where P vi is the intervascular blood pressure (IBP) in (Pa), π vi is the osmotic pressure of the plasma in (Pa), π i is the osmotic pressure of the interstitial fluid in (Pa), L vi is the hydraulic conductivity of the blood vessel walls in (m/Pa s), σ i is the average osmotic reflection coefficient for plasma proteins, and S vi V i is the surface area of the blood vessels per unit volume of tissue in (1/m). The lymphatic drainage term φ li is given as 27,30 where P li is the hydrostatic pressure of intra-lymphatic in (Pa), L l is the hydraulic conductivity of the lymphatic wall in (m/(Pa s)), S li V i is the surface area of the lymphatic vessels per unit volume of tissue in (1/m), and L l S l V is the lymphatic filtration coefficient in (1/(Pa s)). At steady state, the mass continuity equation for incompressible flow in porous media reduces to Now, by combining Darcy's law (1) and the mass continuity Eq. (5), we get where ∇ 2 is the Laplacian operator. The lymphatic drainage term is equal to zero in the tumor region due to the lake of the lymphatic system. Moreover, the mass fluid source term is equal to zero in the ablated-tumor region because the thermal ablation destroys the vascular network. Thus, Eq. (6) can be rewritten as Equation (7), together with the boundary conditions (23)-(34), can be solved analytically or numerically. However, the analytical derivation could be complicated for three regions even there is an analytical solution for www.nature.com/scientificreports/ two regions model, and thus we solve it using COMSOL Multiphysics software. The obtained interstitial fluid pressure and velocity will be used in the drug transport model in the next subsection.
Interstitial drug transport. In this model, the drug source is a miniaturized implant loaded with DOX anticancer drug inserted inside a solid tumor. The implant will release the anticancer drug, which diffuses through the surrounding tissue, and it will be influenced by the various pharmacokinetic processes, including drug efflux/ influx into cells, drug binding/unbinding with interstitial proteins, elimination into blood capillaries, and biodegradation. These factors have a significant impact on drug transport and therapeutic efficacy. Transport of free-DOX in the interstitial space can be mathematically modelled using the following diffusion-convection-reaction equation.
where ∇ is the Del gradient operator and the index i ∈ {a, t, n} refers to the ablated tumor, the non-ablated tumor, and the normal tissues. The function C exi is the spatiotemporal extracellular concentration of free-DOX in the interstitial space in (g/m 3 ), and D Fi is the diffusion coefficient of free-DOX, which has a different value in each region.
The loss rate of free-DOX due to drainage in the lymphatic vessels and elimination by the blood capillaries.
The loss rate of free-DOX due to lymphatic drainage (φ li ) is given by Eq. (4). Since lack of lymphatic system in the tumor region 1,28,30 , we set the loss rate due to lymphatic drainage equal to zero in both the ablated and nonablated tumor regions. Moreover, the experimental studies showed that thermal ablation destroys the vascular structure inside the ablated tumor zone, and therefore the DOX loss into the blood in this region is neglected. However, several blood microvessels may appear in the peripheral region of the ablation zone 37 . Moreover, the initial concentration of DOX in the plasma is assumed to be zero, where the implant is the only drug source 38 .
The loss rate of free-DOX due to elimination by blood vessels (γ vi ) can be expressed as follows 38,39 .
where P Fi is the permeability coefficient of the blood vessel wall to free-DOX in (m/s). The term F B accounts for the binding and unbinding of DOX with Albumin proteins in the interstitial space.
where k a and k d are the DOX-protein association and dissociation reaction rates, respectively, and C bi is the spatiotemporal bound-DOX concertation in the tissues. The doxorubicin molecules can transport to/from the interior of the cell across the cell membrane. Thus, the effect of cellular uptake (influx) and efflux is modelled using the following cellular influx/efflux rate: where D c is the cancer cells density in the unit of (10 5 cells/m 3 ). The cellular influx and efflux functions, ζ inf and ζ eff , are given in Eqs. (14) and (15), respectively.

Modelling of the target cells. Drug absorption by the cells.
The doxorubicin can transport across the cell membrane via passive diffusion and carrier-mediated transport 40 . Free and bound DOX can cross the tumor cell membrane. The amount of bound DOX which enter the cells is neglected 41 . Therefore, we assume that only free-DOX can uptake by the cells, and therefore the intracellular concentration is a function of the extracellular free-DOX concentration. The intracellular DOX concentration C c is expressed in the unit of (ng/10 5 cells), and it changes with the time according to the following equations 42,43 . where ζ inf and ζ eff are the cellular influx and efflux functions, V max is the maximum rate of transmembrane transport, and φ e is the extracellular volume fraction. The parameters V max , k e , and k c were obtained in the literature 43 , by the fitting of the experimental data of intracellular Adriamycin concentration in tumor cells 44 . The parameters mentioned above are listed in Table 2.

Scientific Reports
| (2020) 10:19366 | https://doi.org/10.1038/s41598-020-76123-0 www.nature.com/scientificreports/ Albumin is moving protein molecules that can be found in the bloodstream and interstitial space. Some free-DOX molecules may bind to proteins, such as Albumin, in the interstitial space, then new macromolecules (bound-DOX) will be created. The bound-DOX may unbind from the DOX-protein complexes and then becomes free. The spatiotemporal transport of bound-DOX in the extracellular space can be modelled using the following diffusion-convection-reaction equation 27,41 : where C exi is the spatiotemporal extracellular concentration of bound-DOX in the interstitial space and D Bi is the diffusion coefficient of bound-DOX, which has a different value in each region.
The loss rate of bound-DOX due to the blood and lymphatic microvessels is given as where P Bi is the permeability coefficient of the blood vessel wall to bound-DOX.
Death of the tumor cells. The cancer cell density can be described using the following equations, which include the impact of the free-DOX concentration and the natural growth and death of cells 27,42 .
where k gr and k dr are the natural growth and decay rate constants of the tumor cells, respectively, and k m is the saturation constant. The nonlinear function K reflects the effect of the anticancer drug, which depends on the intracellular concentration of free-DOX, the maximal DOX cell killing rate ( k max ), and Michaelis constant EC 50 . Initial tumor cell density is used as an initial condition for solving Eq. (18). The parameters mentioned above are listed in Table 2.
Drug implant model. In this study, the drug implant is loaded with anticancer drug DOX. The main design parameters of the implant are the release rate, the implant size, and the amount of loaded drug. The release rate depends on the implant formation and the physicochemical properties of the loaded drug, which can be adjusted during the design phase by selecting appropriate materials, e.g., polymer and drugs 48 . The optimal release rate can help in improving the therapeutic outcomes by reducing the side effects, which in order save time and cost. The amount of the released drug can be experimentally monitored over time to obtain the release profiles. Relevant mathematical models for the release kinetic can be fitted to the experimentally measured release curves to predict the release rate constants 49 . The implant has a dual-release pattern, i.e., it releases DOX over two phases: fast burst release over a short time duration followed by slow sustained release over an extended period of time.
Thus, we can model both dual-release and sustained release implants by adjusting the release parameters. For example, in-situ forming implants (ISFIs) and double-layer implants have dual-release pattern 6,12 .
The cumulative amount of released DOX at the time t can be mathematically modelled using the bi-exponential first-order kinetic model as www.nature.com/scientificreports/ where M 0 is the total amount of loaded-drug in the implant in (mg), W ∞ is the total fraction of drug released at steady state, and f is a fraction of drug released during the burst phase. The parameters k f and k s are release rate constants for burst and sustained release phases in (s −1 ), respectively. Now, the rate of drug release across the spherical implant surface at time t per unit area of the implant surface in g/(s m 2 ), i.e., the flux, can be expressed as where A = 4πr 2 0 is the surface area of the implant and r 0 is the implant radius. The characteristic parameters of the dual-release implant are listed in Table 3. Similar to the experimental release data 6 , we chose 10% of the loaded drug to be released within the first hour while different values of sustained-release rate constant are examined.
The experimental studies in the literature showed insignificant variation for the size of the implant during the release duration since the implant releases drug prior to any significant degradation 38 . Thus, it is reasonable to assume that the size variation due to biodegradation is negligible during the observation time 6,38 . Model parametrization. In this work, the values of the model parameters are obtained from the published experimental data and other studies in the literature, assuming that the growth of the tissues is negligible during the observation timeframe. This assumption is widely used in the literature since the time scales of fluid and drug transport phenomena are relatively short 27 . Thus, it would be reasonable to assume that the system's physiological parameters to be time-independent. The parameters are defined through the paper and summarized in Tables 2, 3, 4 and 5.
Microvasculature density (S/V). The microvasculature density is the ratio of vascular surface area per unit volume of tissue. The microvasculature density of capillaries highly varies among tumor types and within the same type of tumor 34 . However, the experimental data show a leak of lymph vessels in the tumor tissues. Different studies show that the larger tumors have smaller microvasculature densities 34,50 . For a tumor with 2 cm diameter (i.e., volume = 4188 mm 3 ), the surface area of the blood vessels per unit volume of tumor tissue is approximately equal to 10 4 1/m 34 . In normal tissues, the surface area of the blood vessels per unit volume of tissue is measured as 7 × 10 3 1/m 51 . As mentioned before, the histological analysis confirms that the thermal ablation destroys the vascular network and tumor cells in ablated tumor tissue 4,8,12,21,[52][53][54] . Thus, the effect of microvasculature density and, consequently, drug loss through blood microvessels in ablated tumor tissue is neglected. www.nature.com/scientificreports/ Extracellular space fraction (φ e ). The volume fraction of extracellular space in the tumor is much larger than that in the normal tissue, and it ranges from 0.2 to 0.6 with tumor cell density ranges from 0.955-15.3 10 5 cells/ mm 341 . In this study, the volume fraction of extracellular space φ e and the initial cell density D c0 are chosen to be equal to 0.4 and 10 × 10 5 cells/mm 3 , respectively.
Microvasculature permeability (P F , P B ). The vasculature permeability measures the capability of the blood or lymph microvessels to exchange various substances in and out of the vasculature. The vasculature permeability of Albumin (similar to Albumin bound-DOX) and free-DOX is approximately threefold higher in tumor than that in normal tissues 41 .

Diffusion coefficients (D F , D B ).
The diffusivity of free-DOX is higher than Albumin-bound DOX, where the molecular weights (MW) of free-DOX and bound-DOX are 544 Da and 69 kDa, respectively. Moreover, the diffusivity of free and bound-DOX in the tumor is larger than that in normal tissues, as listed in Table 5. However, the experimental studies found that the diffusivity of DOX near the center of the ablated liver tumor after RFA, D ac , is 75% higher than that in non-ablated tumor 20,23 . Based on the histological findings in ablated tumor tissues, the diffusivity shows dependency on the radial distance with a higher value at the center of the ablated tumor than the periphery region. The diffusion coefficient in the outer region of the ablated tumor ( r c ≤ r ≤ r a ) is characterized as 20,23 where r a is the ablation zone radius, r c = αr a is the radius of the central region of the tumor, and α = 0.47.
The diffusivity D ac of the tumor center decreases linearly with the radial distance to finally reaches the diffusivity of the non-ablated (risk) tumor, i.e., D t . In this model, we use an average diffusion coefficient within the ablated liver tumor tissue. The mean diffusivity of the outer ablated region is calculated using the scale relationship of DOX diffusivity in ablated and non-ablated tissues, then by taking the average of the diffusivity given by Eq. (22). Model geometry. In this work, a 3-D spherical solid tumor is considered with a diameter of 2 cm, as shown in Fig. 11. This value falls within the range of different tumor sizes encountered in reality 59 , e.g., the diameter of tumors in rats and rabbits ranges from 0.5 to 2 cm. We examine three cases for the tumor microenvironment, namely, non-ablated tumor, 80% ablated tumor, and 90% ablated tumor. For example, 80% ablated tumor means that a tumor region of a radius ( 0.8 × r t ) is ablated. In the case of the ablated tumor, a thin rim of viable cancer cells (risk region) will remain at the tumor periphery. The thickness of the risk regions with 80% and 90% ablation zones are 2 mm and 1 mm, respectively. Here, the thickness of the normal region is chosen to be 5 mm. The implant has a radius of 1.5 mm, and it is located at the tumor center. The model geometry and mesh are created using a built-in CAD kernel in COMSOL Multiphysics package. The mesh size is chosen as "Finer" based on a convergence mesh independence test, which shows that a 5-times decrease in the mesh element size will provide a negligible enhancement, i.e., < 5%, in DOX concentration profiles. Moreover, we refined the mesh using "Boundary Layer Setting" at the implant/tissue boundary and other interface boundaries between the various regions, to handle the rapid change of drug concentration, velocity, and pressure at these interface layers.
Boundary conditions. Due to spherical symmetry, the pressure at the tumor center is characterized using no-flux boundary condition as In this work, the observation time scale in the numerical analysis is assumed to be much shorter than the time scale for the growth of the tumor and normal tissues 27 . Therefore, the interface boundary between the various regions are assumed to be fixed. The continuity boundary conditions of the interstitial pressure and fluid flux are imposed at these interface boundaries as  www.nature.com/scientificreports/ In addition, the outer boundary of the normal region is assumed to be fixed. The interstitial fluid pressure at this boundary has a minimal constant value. Thus, Dirichlet boundary condition can be applied at the outer boundary as The total interstitial drug flux is a combination of diffusion and convection fluxes. In this model, the interstitial velocity field mainly appears at the interface boundary between the tumor and the normal tissues due to a large pressure difference at that layer. The continuity boundary conditions of drug flux and concentration are applied at the interface boundaries between the various regions. The continuity boundary conditions are applied for both free and bound DOX as where ∇ is the gradient operator in the spherical coordinate system and the index j is used to indicate either the free or bound DOX.
No flux (Neumann) boundary condition is applied at the outer boundary of the normal region as The release process of DOX drug from the implant surface is modelled using the following flux boundary condition: where F rs (t) is given by Eq. (21).