Spatiotemporal multi-scale modeling of radiopharmaceutical distributions in vascularized solid tumors

We present comprehensive mathematical modeling of radiopharmaceutical spatiotemporal distributions within vascularized solid tumors. The novelty of the presented model is at mathematical level. From the mathematical viewpoint, we provide a general modeling framework for the process of radiopharmaceutical distribution in the tumor microenvironment to enable an analysis of the effect of various tumor-related parameters on the distribution of different radiopharmaceuticals. We argue that partial differential equations (PDEs), beyond conventional methods, including ODE-based kinetic compartment modeling, can be used to evaluate radiopharmaceutical distribution in both time and space. In addition, we consider the spatially-variable dynamic structure of tumor microvascular networks to simulate blood flow distribution. To examine the robustness of the model, the effects of microvessel density (MVD) and tumor size, as two important factors in tumor prognosis, on the radiopharmaceutical distribution within the tumor are investigated over time (in the present work, we focus on the radiopharmaceutical [18F]FDG, yet the framework is broadly applicable to radiopharmaceuticals). Results demonstrate that the maximum total uptake of [18F]FDG at all time frames occurs in the tumor area due to the high capillary permeability and lack of a functional lymphatic system. As the MVD of networks increases, the mean total uptake in the tumor is also enhanced, where the rate of diffusion from vessel to tissue has the highest contribution and the rate of convection transport has the lowest contribution. The results of this study can be used to better investigate various phenomena and bridge a gap among cancer biology, mathematical oncology, medical physics, and radiology.

Tumor heterogeneity is among the most important contributing factors to cancer treatment failure [1][2][3][4][5] . Precision oncology aims to address this issue by investigating tumor biology variations and tumor micro-environment (TME) heterogeneity. Personalized treatment planning is required to mitigate the effect of these variations, away from a one-fits-all approach. The theranostic paradigm is one of the best platforms to actualize precision oncology [6][7][8][9] . Its therapeutic aspect, through radiopharmaceutical therapy, targets specific tumor cells while its diagnostic counterpart helps to illustrate tumor heterogeneity. This combination is a perfect match for personalized cancer care. In this paradigm, similar radiopharmaceuticals are used for diagnosis and therapy. Radiopharmaceuticals are labeled with two categories of radionuclides, one for diagnostic purposes and the other for treatment. The diagnostic radionuclide usually is a positron emitter (such as 18 F, 13 N, 11 C, or 15 O), while the therapeutic radionuclide is an alpha or beta emitter (such as 225 Ac, 177 Lu, or 90 Y). The critical property of these two radiopharmaceuticals (i.e., diagnostic and therapeutic) is their pharmacokinetic similarities. This phenomenon will provide physicians with a unique opportunity to image the tumor heterogeneity on the one hand and to address this depicted heterogeneity through therapy adjustment on the other. Tumor heterogeneity has two components: tumor cell variations (inter-tumoral and intra-tumoral mutations) and tumor micro-environment heterogeneity (including vascularity, permeability and diffusability, hypoxia, and the extend of immune system presence and function, among the others) [10][11][12] . Quantitative imaging using dynamic positron emission tomography (PET) represents the aggregation of these two factors. Traditionally, the second category (tumor micro-environment) was assumed insignificant, and the quantitative models did not consider spatial heterogeneity of the micro-environment. This assumption appears unrealistic in light of recent discoveries about the pivotal role of the tumor micro-environment in treatment success and failure 13,14 . To address this unmet and critical need, the conventional temporal models (such as compartmental kinetic analysis based on ordinary differential equations (ODEs) 15,16 ) should be updated using spatiotemporal distribution models (SDMs) based on partial differential equations (PDEs)) [17][18][19][20][21][22][23][24] . This study aimed to model one of the microenvironment variables: neovasculature.
The majority of therapeutic radiopharmaceuticals in current medical practice are small molecules, in contrast to large molecules such as monoclonal antibodies (mAb). Since the pharmacokinetics of these two groups are significantly different, we decided to select a small molecule. Considering the availability of experimental data, clarity of biochemical processes, and relevance to tumor biology, [ 18 F]-Fluorodeoxyglucose ([ 18 F] FDG) was selected to represent small molecule kinetics, while our intention is to extend and apply this framework to theranostics in future efforts.
Several decades ago, biochemical studies in tissue cultures demonstrated that metabolic abnormalities, especially malignant tumors, show high rates of glucose uptake and glycolysis 25 . [ 18 F]-Fluorodeoxyglucose ([ 18 F] FDG) PET allows the use of these metabolic abnormalities of tumors for clinical diagnosis 26,27 . The FDG uptake into tissue reflects the transport and phosphorylation of glucose by viable cells, where FDG is transported across the cell membrane and subsequently is phosphorylated by hexokinase. In contrast to glucose, FDG-6-phosphate cannot be further metabolized; moreover, because FDG-6-phosphate is a highly polar molecule, it cannot diffuse out of the cell and remains intracellular. Therefore, due to this trapping mechanism, the concentration of FDG is steadily increasing in metabolically active cells, and as a result, the contrast between the tumor and normal tissue increases.
The angiogenesis process, formation of new blood vessels from pre-existing blood vasculature, is a crucial component for supporting tumor growth by providing a source of nutrients and oxygen 28,29 . When the volume of the primary tumor exceeds a threshold ( ∼ 1 mm 3 ) 30 , its central cells become deficient in oxygen (hypoxic state); in response, the tumor and the cells surrounding secrete a variety of chemotactic and morphogenic tumor angiogenesis factors (TAFs) such as vascular endothelial growth factor (VEGF) 31 . The binding of TAFs to its receptor on the ECs of the existing vessel wall leads to the differentiation of ECs and the formation of tip endothelial cells (tECs) with increased motility. The tECs migrate toward the tumor in response to the TAF gradient (i.e., chemotaxis) and adhesion gradient to fibronectin (i.e., haptotaxis) 32 . Extend of tECs from the existing vessel wall, create new blood vessel sprouts in the initial stages of angiogenesis 33 . Transverse movement of the tECs leads migrating vessels to meet and forms connections to create a complex network of immature microvessels that supports blood flow into the tumor. Therefore, the process of tumor-induced sprouting angiogenesis is governed by the coordinated cell dynamics, which involves initial sprouting, branching, anastomosis, and vessel remodeling 34 .
Current mathematical simulation models of tumor angiogenesis can be classified into three categories: continuous, discrete, and hybrid models. The continuous approach assumes the vascular cells as bulk and continuous aggregates without considering the cells alone [35][36][37][38] . Although this model provides good qualitative results, it can not model the morphological details of the vascular network structure. The discrete models assume a set of defined rules and building-blocks structure very close to reality to simulate EC migration and the creation of new vessels [39][40][41] . Hybrid models consist of a combination of continuous and discrete models that trace the EC pathways using variable motion probabilities [42][43][44][45] .
In this study, we present the design and implementation of a comprehensive spatio-temporal multi-scale framework to collectively model the radiopharmaceutical distribution for different tumor sizes and microvessel densities (MVD). Our model includes subcellular, cellular, and tissue level size scales. The subcellular scale consists of biochemical agents fibronectin and tumor angiogenesis factor. The cellular scale includes the movement of endothelial cells and uptake of radiopharmaceuticals into the tumoral and normal cells. At the tissue scale, www.nature.com/scientificreports/ blood vessel growth is implemented by a hybrid method for the movement of tECs in two-dimensional interstitial space. We also model dynamic simulation of blood flow/pressure, which is determined by mechanobiological and biochemical signals from wall shear stress, pressure, and metabolic stimuli with accurate hemorheology and hemodynamics at the tissue level. In addition, the interstitial fluid flow that influences the distribution and transport of radiopharmaceuticals is also calculated. Subsequently, radiopharmaceuticals transport in different scales is analyzed using SDM equations.

Methods
For our modeling, we assume a semi-realistic micro-scale architecture of vasculature, interstitium, and tumor, as schematically depicted in Fig. 1. In this structural schema, various biochemical (e.g., TAF and O 2 ) have spatially heterogeneous concentrations, and their gradients are highlighted at the border. In the following section, first, we present a mathematical model of angiogenesis and dynamic simulation of blood flow. Subsequently, spatiotemporal modeling of radiopharmaceutical transport in different scales will be described. Computational implementation will be discussed at the end (including geometry, computational domain, boundary conditions, and other simulation parameters).
Mathematical model of angiogenesis. The present angiogenesis model is inspired by the tumorinduced angiogenesis model initially proposed by Anderson and Chaplin 43,45 . The mathematical model predicts capillary formation by tracking the motion of ECs at the tip of the capillary sprout, ultimately forming a vascular network. The details of the rules have been determined in the previous publication 43, 45 . In summary, three main mechanisms for the motion of ECs are considered in this model: random motility, chemotaxis, and haptotaxis. Each motion component is the function of the motion stimulus gradient. Random motility is modeled like diffusion of mass due to the concentration gradient. Chemotaxis is related to the TAF concentration gradient. TAF diffuses through the interstitium and activates the ECs. The distribution of TAF has its maximum value near the tumor zone and decreases to its minimum near the parent vessel (Fig. S1). Haptotaxis is the motion of cells due to the adhesion gradient in the environment. Cells move in search of stronger bonds in the interstitium, which is expressed in relation to the fibronectin concentration gradient in the environment. The complete system of equations describing the interactions between EC density, n , fibronectin, f and TAF, c TAF , are respectively written in Eqs.  is determined in five states by P i coefficients, which include the probability of the cell being stationary ( P 0 ), or moving right ( P 1 ), left ( P 2 ), up ( P 3 ), or down ( P 4 ). Anastomosis and branching occur during the angiogenesis process. It should be mentioned that Microsoft Office PowerPoint 365 was used to create this figure.
D n is the random motility coefficient of EC, χ is the chemotaxis coefficient, α is constant of chemotaxis coefficient, φ is haptotaxis coefficient, β and γ are production and uptake coefficient of fibronectin, respectively, and η is uptake coefficient of TAF by ECs. Implementation of the angiogenesis method to calculate the movement probabilities for each endothelial cell derived from these equations is presented in Supplementary Material (Eqs. (S1-S8)).

Hemodynamics and interstitial fluid flow.
For the fluid dynamics simulation model, we consider flow in three regions: intravascular blood flow, transvascular fluid flow, and interstitial fluid flow. To study the incompressible flow through the microvascular network, the flow rate in each vessel is calculated by applying mass conservation at each network junction. The continuity equation for intravascular blood flow at an interconnecting point, like c, in the network is written as Eq. (4) (Fig. S3). In this equation, the index k refers to adjacent nodes, Q k c is the net blood flow rate for each capillary calculated as the difference between the intravascular blood flow rate, Q k b.c , and the transvascular fluid flow rate, Q k t.c . The Womersley flow number for blood flow in the capillary network is very low due to the very small diameter of the newly created vessels. Hence, the pulsatile effects of the cardiac cycle can be ignored. Thus, Hagen-Poissville's law can be applied for intravascular blood flow as an exact solution to the Navier-Stokes fluid dynamics equation (Eq. 5) 41 . Transvascular fluid flow rate is calculated by Starling's law, which indicates the role of oncotic and hydrostatic pressures in fluid movement across capillary membranes (Eq. 6) 41 .
In Eq. (4), N is the number of peripheral vessel lattice nodes adjacent to the central vessel node. In the 2-D simulation for a fully connected network, N is 4 and β k c is a positive integer '0' or '1' , which describes whether nodes k and c are connected ( β k c = 1) or not ( β k c = 0). In Eq. (5), P b is intravascular pressure (IVP), H is hematocrit, and d and l are the diameter and length of the new vessel, respectively. In Eq. (6), L p is the hydraulic conductivity of the microvascular wall, P b is the average IVP in each element, P i is the average interstitial fluid pressure (IFP) outside of the vascular element, π b is the osmotic pressure of the intravascular plasma, π i is the osmotic pressure of the interstitial fluid, and σ is the average osmotic reflection coefficient for plasma proteins.
IVP is calculated by applying Poiseuille's equation for the flow in the vessels. In contrast, interstitial pressure for the peripheral tissue of a microvascular network is calculated by solving the governing equation for fluid flow through a porous medium. The mass balance equation for a steady-state incompressible fluid is modified by adding source and sink terms for biological tissues (Eq. 7). Blood vessels are fluid source terms, and lymphatic vessels are sink terms, which are respectively shown in Eqs. (8) and (9), according to Starling's law. Moreover, Darcy's experimental observations show that the interstitial fluid velocity (IFV) in porous media is proportional to the pressure gradient (Eq. 10).
v i is the IFV, φ b is the rate of fluid flow per unit volume from blood vessels into the interstitial space, φ L is the rate of fluid flow per unit volume from the interstitial space into lymph vessels, S/V is the surface area of the (1) ∂n ∂t = D n ∇ 2 n random motility − ∇ · χ(1 + αc TAF )n∇c TAF + φn∇f chemotaxis and haptotaxis www.nature.com/scientificreports/ microvascular per unit volume for mass transport in the interstitium, P L the hydrostatic pressure of the lymphatic, and K is the hydraulic conductivity of the interstitium. By combining Eqs. (7) and (10), we can derive the Poisson-Laplace's equation to calculate IFP in both normal and tumor tissues (Eq. 11).
Hemorheology. Blood has significant non-Newtonian properties in low Reynolds numbers. We assume that blood is a suspension of non-Newtonian red blood cells within a Newtonian plasma fluid [46][47][48] . Pries et al. 49 presented an empirical relation for the dynamic blood viscosity in capillaries, shown in Eq. (12). In Eq. (12), µ blood , µ plasma and µ rel are the dynamic viscosity of blood and plasma and apparent viscosity of blood, respectively. The apparent blood viscosity is defined as a function of vessel diameter and hematocrit by Eq. (13). µ 0.45 is the relative apparent blood viscosity for a fixed hematocrit of 0.45, which is defined in Eq. (14) as a function of vessel diameter, and C is a function that describes the shape of viscosity dependency on hematocrit, which is defined in Eq. (15). In Eqs. (13)(14)(15), the unit of d is the µm.
Hematocrit distribution has a significant role in simulating the hemodynamic characteristic of the microvascular network, which depends on blood flow characteristics such as velocity. In general, the hematocrit distribution at vessel bifurcations can change depending on the flow velocity in each branch. In other words, if the velocity ratio of the two branches exceeds the threshold value, U cr , all hematocrit enter the faster branch at bifurcations 50 . According to this, the relation between the hematocrit of the parent vessel, H i , and the branches, H 1 and H 2 , are written as Eqs. (16)(17)(18), respectively, based on the velocity ratio of the two branches, U 1 /U 2 . is a phenomenological parameter that accounts for the strength of the non-symmetry of the hematocrit distribution at bifurcations 50 .
Dynamic structure adaptation method. Capillaries are able to continuously adapt their diameter in response to physical and biochemical stimuli of the tissues that the capillaries supply 47,49 . The wall shear stress stimuli, S wss , and the IVP created by blood flow,S p , and the biochemical stimuli such as the hematocrit-induced metabolic, S m , lead to the remodeling of the vascular diameter, shown in Eqs. (19)(20)(21), respectively. τ w is the wall shear stress in a capillary vessel defined by Eq. (22). Experimental observations show that when each vessel senses the shear stress of the wall, it adjusts its diameter to achieve a uniform level of stress. τ ref is a small constant to avoid singular behavior at low wall shear stress rates. τ e (P b ) is the wall shear stress resulting from the blood pressure, P b , defined in Eq. (23). Q b is the rate of blood flow in the vessel, and Q ref is the maximum value of Q b within the network.
Normal tissue Tumor tissue For existence blood source Normal tissue Tumor tissue Otherwise www.nature.com/scientificreports/ Finally, the change of diameter for each vessel in the network in proportion to the total signal and time step is considered in Eq. (24). k p is the adaptive response sensitivity of the vessel diameter to changes in IVP, k m is the adaptive response sensitivity of the vessel diameter to changes in metabolic state, and k s is the shrinking tendency of the vessel in the absence of positive growth stimuli.
Spatiotemporal distribution of FDG. In PET, images are a composite of different superimposed signals. A signal, for example, may describe the amount of the phosphorylated radiopharmaceutical. In order to separate the desired component of the signal, a mathematical model is required that relates the dynamics of all radiopharmaceutical states to the resulting PET image. Mathematical kinetic models are used to analyze the time sequences of PET images, in which each of the radiopharmaceutical states is known as a compartment. Each compartment is determined by the concentration of the radiopharmaceutical inside it as a function of time. These concentrations are related through a set of ODEs, which express the balance between the input and output of each compartment. In compartmental modeling, it is assumed that there are no spatial concentration gradients within the sampled area (e.g., a voxel). The directly measured radiopharmaceutical concentration in blood as a function of time acts as a model input function. The coefficients of the differential equations in the model are considered as constants that reflect the inherent kinetic properties of the specific radiopharmaceutical molecule in the system. By formally comparing the model output to the experimentally obtained PET data, the values of these kinetic parameters can be estimated, and information about the delivery process can be extracted 51 . The 4-compartment 5-rate-constant (5 K) model distinguishes the kinetic steps of [ 18 F]FDG delivery to the extracellular matrix (ECM), its transport from the extracellular to the intracellular space, and its intracellular phosphorylation 16 . The 5 K model is described by ODEs, shown in Eqs. (S9-S11) at the Supplementary Material, and can be viewed as an extension of the classic model by Sokoloff et al. 15 . In fact, a value of the 5 K model is in its explicit accounting of an extracellular compartment.K 1 and K 2 are exchange rate parameters between plasma and extracellular space, K 3 and K 4 are transport rate parameters into and out of the cell, and K 5 is the phosphorylation rate.
As mentioned before, in contrast to classic compartment model, the SDM model uses PDEs that enables the evaluation of radiopharmaceutical distribution in both space and time scales, hence distribution in space is not modeled as independent from time. Fig. 2 describes the SDM model. From the arterial blood, the radiopharmaceutical passes into the second compartment, known as the extracellular compartment, where it is considered that the radiopharmaceutical is redistributed in both time and space. The third and fourth compartments are the region of intracellular radiopharmaceutical. Movement between the plasma, extracellular, and intracellular spaces is governed through five parameters: L 1 , L 2 , L 3 , L 4 , and L 5 , similar to the rate constants K 1 , K 2 , K 3 , K 4 , and K 5 seen in conventional compartmental models. The general form of the SDM equations in tissue is shown by Eqs. (25)(26)(27) as 52 :  , v i is the IFV vector,σ f is the filtration reflection coefficient, P m is the vascular permeability coefficient, and Pe represents the Peclet number, indicating the ratio of convection transport rate to diffusion transport rate across the vessel wall to interstitium. The first and second terms in the right-hand side of Eq. (25) determine the diffusion transport and convection transport of radiopharmaceuticals within the interstitium, respectively.

Computational implementation. Computational geometry and simulation cases.
To reduce computational cost and simplify the calculations, we selected a L × L square as the computational domain representing TME and its surrounding normal tissues, as demonstrated in Fig. 3. The "D" in this figure shows the tumor size and is non-dimensionalized by the domain length ( D L = 0.2 ). A parent vessel is located the left edge of the domain where the vascular network starts to grow. To investigate the effect of tumor size, three different tumor sizes including D, 2 × D , and 3 × D are considered. On the other hand, to examine the impact of MVD within the TME, two different states are considered for initiating angiogenesis from the parent vessel, including network with 3 initial sprouts and network with 5 initial sprouts, at the same distance from each other. The values of the model parameters used in this study are listed in Table S2 of the Supplementary Material.  www.nature.com/scientificreports/ Solution strategy. In this study, a two-dimensional hybrid sprouting angiogenesis method is applied to produce the microvascular network. The governing equations are discretized by a finite difference scheme to a set of algebraic equations. The solution of the discretized equations is performed in an iterative process by MATLAB. The calculation of blood flow in capillaries includes a set of non-linear equations. Hence, an iterative method is applied to solve the blood flow equations in the capillary network by MATLAB (160,000 elements). Subsequently, the momentum and mass equations in the microvascular network and interstitial space are solved. The resulting IFP and IFV values are used for solving the SDM equations to obtain the distribution of different concentrations (28,082 elements). The initial and boundary conditions of this study are mentioned in Supplementary Material (Table S1). A step-by-step flowchart of the present multi-scale model is reflected in Fig. 4. The Darcy, and SDM equations are solved by the commercial finite element software COMSOL Multiphysics version 5.5a (COMSOL Inc., Stockholm, Sweden). In addition, the residual square errors are set to 4 orders of magnitudes and the time of analysis for the FDG distribution is considered 1 h. All simulations described in this manuscript were performed on an Intel Core i5-8250U, 1.8 GHz CPU and RAM 8 GB computer.

Results and discussion
Morphology of the angiogenesis network. To analyze the different stages of tumor angiogenesis, a sample vascular network induced by two parent vessels on the left and right sides of the domain has been generated over 30 days (Fig. 5). At 5 days, the TAF forms a gradient in the domain that reaches the parent vessels and causes sprouting angiogenesis to form new vessels. Between days 5-15, the generated vessels are elongate, branch, deform and spread toward the tumor. After day 15, new vessels have been extended within the TME to reach together and form vascular loops. A qualitative comparison between the results of this study and biological observations of tumor angiogenesis morphology 53,54 is shown in Fig. 6. In agreement with in vivo observations, the simulations show that the branching of new vessels increases near the tumor boundary as well as within the tumor where the TAF concentration is high (Fig. 6a-c), because TAF induces new tECs and create new branches. As mentioned, to reduce computational  www.nature.com/scientificreports/ cost and simplify the calculations, the angiogenesis model caused by one parent vessel has been considered in this study (Fig. 6d). This assumption is biologically validated and confirmed by Fig. 6e,f. In the calculations, given the sake of investigating the dependency of FDG radiopharmaceutical activity on the tumor size as well as MVD within the tumor, six capillary networks-including three different tumor sizes (D, 2 × D , and 3 × D ) in two different states of initiating angiogenesis (with 3 initial sprouts and 5 initial  www.nature.com/scientificreports/ sprouts)-have been generated and examined in detail, as demonstrated in Fig. S2 of Supplementary Material. It should be noted that in a realistic capillary network, the maximum distance between the tumor cells and the closest microvessels is considered to be about 100 μm 55 . This is not considered in our microvascular networks, to more clearly investigate the effect of MVD and the distribution of solute within the tumor. However, this assumption satisfies the mechanisms of extravasation from microvessels and transport in interstitium. In fact, MVD plays a potential role as prognostic biomarker for different cancers including breast 56,57 and has been shown as elevated in higher-grade tumor tissues compared to lower-grade ones 58 64 , which demonstrated the range of 586-4200 Pa for IFP in the tumor. Moreover, the IFP value in normal tissue is around 53.7-97.5 Pa, that agrees with the values of the experimental study 65,66 , which demonstrated the range of − 400 to 800 Pa for IFP in normal tissue. IFP has its greatest value in the tumor region because the lymphatic system is dysfunctional and the capillary network of the tumor has a higher leakage rate than that of normal tissue, leading to radiopharmaceutical accumulation in the tumor interstitium, and as a result, an elevated IFP.
IFP is known as an important barrier against efficient solute transport to tumor tissue. IFP is proportional to the MVD in tumor and normal tissues. Hence, the IFP value is higher in the area where the capillaries are closer together. In addition, tumor with higher MVD has greater IFP, according to the mean values of the measured IFP. As Fig. 7 illustrates, the value of IFP decreased rapidly at the tumor boundary because the transvascular flow is absorbed by the lymphatic vessels in the normal tissue region. Moreover, the heterogeneous capillary network, as the source term in the interstitial fluid flow equation, causes the heterogeneous distribution of IFP in tumor tissue.
Given a direct relationship between IFP and IFV based on Darcy's law, IFV distribution can be obtained in the whole tissue domain, as demonstrated in Fig. 8 for six considered networks. The velocity has very low values in most areas. The fluid flow direction is from high pressure area toward a low pressure area. Indeed, the maximum   67 , and experimental observation of Hompland et al. 68 . In the following, the effect of five factors on flow simulation is examined: Osmotic pressure. The transvascular flow, which is of crucial importance for the homeostasis of tissues, is determined by the hydrostatic and osmotic pressures in vessels and in the interstitium. In tumors, normal osmotic pressures necessary for the transit of solutes into and out of vessels along gradients are severely compromised owing to the local tumor microenvironment 69 . Hence, the hydrostatic pressure difference between intra-and extravascular spaces is the dominant term of transvascular flow 70 . Therefore, due to the minor role of osmotic pressure in the transvascular flow, it can be neglected 71 , although this assumption is far from the truth.
Hematocrit distribution. As can be seen from the equations, intravascular flow ( Q b ), dynamic blood viscosity ( µ blood ), and metabolic stimuli for changing the diameter ( S m ) of each branch depend on the hematocrit in that branch. No change in hematocrit of the branches affects the mentioned cases. According to the researches on vascular dynamics, the hematocrit distribution in branches varies proportional to velocity.
Changing vessel diameter. If we call a network with change in diameter a "dynamic network" and a network with fixed diameter a "static network", the dynamic network adapts its vessels diameter to reach enough blood flow in each vessel. Therefore, the intravascular pressure of the dynamic network is higher than the static network (Fig. S5). In addition, the dynamic network has a higher IFP because of the elevated intravascular blood pressure in the tumor, which increases the transvascular flow rate. The results of the blood pressure distribution for the static network led to an unrealistic distribution of solute concentrations.
Lymph pressure. Lymphatic vessels are absent or non-functional in tumors 72,73 . According to the equations, if the lymphatic term is considered in tumor, it will act as a sink term. Therefore, it removes both the drug and the fluid inside the tumor, which reduces the IFP and radiopharmaceutical accumulation. Moreover, if P L > 0 , the difference between P i and P L becomes larger and the sink term effect decreases in normal tissue.
The removal of blind ends and void of flow vessels. Since a mathematical method is used to generate the network, all the generated branches may not be placed in a flow loop. These branches, which are generated due to the nature of the mathematical model, are not practical from a physical viewpoint. As stated by 74 and applied in 41,75-77 , elimination of vessels which are void of flow has proven to be applicable in mathematical models. Therefore, the capillaries that have a flow rate of less than 1% of the network's maximum flow rate can be elimi- Distributions of FDG uptake. A qualitative comparison between the results of this study and the real FDG PET image of a large polyp-shaped gastric cancer for a 70-year-old man 80 is shown in Fig. 9. In agreement with the real PET image, the simulations show high uptake of [ 18 F]FDG in the tumor region. In addition, the created color gradient clearly demonstrates the difference between normal and tumor tissue. We emphasize that this comparison with experimental data is only for demonstration purpose and the simulated FDG image (Fig. 9a) does not use any clinical data.
Since the plasma concentration of the radiopharmaceutical varies by time, the arterial concentration curve of FDG, obtained from Backes et al. 81 , is used for the C P concentration profile (Fig. 13). We assumed that C P is uniformly distributed throughout the vascular network and will enter into the interstitium 67 . The spatiotemporal distribution of FDG concentrations at seven different post-injection time frames for six considered networks are shown in Fig. 10 and Figs. S6-S8 of Supplementary Material. At the beginning of FDG infusion, C i is dominant compared to C e and C m . Over time, the value of C i decreases and first C e and then C m dominate the total concentration ( C tot = C i + C e + C m ). Since the vessels in the capillary network act as source terms for the radiopharmaceutical, the maximum concentration values can be found in the vicinity of the capillary walls. The radiopharmaceutical is then transported and distributed in the other parts of the domain via diffusion and convection mechanisms.
According to the results, the networks induced from the 5 initial sprouts show a more uniform distribution due to the higher MVD. For networks induced from the 3 initial sprouts, due to the heterogeneous structure of the capillary network in tumor region, the radiopharmaceutical delivery to some cells is limited, and as such, uptake is reduced. Furthermore, the maximum total concentration occurs in the tumor region at all time frames due to the high permeability of tumor capillaries compared to normal tissue and also the higher rate of metabolism in cancerous tissue than normal tissue. We also considered the 'uniform network' model by assuming uniform solute transport throughout the tumor. In the uniform network model, the vessels do not exist physically. Instead, it is assumed that they are existed in the whole domain as a source term including tumor and normal tissue. The capillaries' effects have been included in the mathematical modeling by S/V parameter in the equations as a source term in whole the domain (in Eq. S12-1 for concentration transport source term and in Eq. (8) for fluid transport source term); therefore, geometry of parent vessel is not considered in this case. Similarly, C P as the injected concentration term, acts as a parameter in the source term equation (in Eq. S12-1 for concentration transport source term) 61,82 . This model can show the effect of the extreme case of MVD on FDG concentration. Figure 11 shows the mean FDG concentration across the tumor tissue of the considered networks over time. The results include the comparison of outcomes of all six investigated cases with the uniform network case. In addition, the results of FDG distribution were averaged spatially and reported.
Equation (25) in the SDM equations consists of three processes related to diffusion transport ,convection transport, and chemical reactions. The results clearly show that MVD of the capillary network has a significant effect on FDG uptake. The uniform network shows a higher concentration value than the vascularized models. In this case, the solute distribution is only dependent on diffusion from vessels due to the uniform distribution of the capillary network; therefore, the effect of transport in tissue is greatly reduced. At all times, the uptake of the dense network induced from the 5 initial sprouts is greater than that of the 3 initial sprouts, at equal tumor www.nature.com/scientificreports/ size. By increasing the MVD of networks, the effect of terms related to spatiotemporal transport decreases and the contribution of diffusion terms found a significant effect. When MVD is high, the contribution of diffusion from vessel in the solute transport is greater than diffusion in tissue interstitum; however, with decreasing MVD, the effect of diffusion in tissue interstitum increases. In other words, the diffusion term from vessels depends on the structure of the capillary network. Additionally, as shown in Fig. 12a,b, as the size decreased in tumors by the same MVD, the mean total FDG concentration increased; the contribution of diffusion from vessels increases by decreases in tumor size. Moreover, as shown in the previous section, in the smaller tumor size, the IFP values are lower, leading to an increase in the convection rate of the source term in Eq. (25). Indeed, by reducing the IFP value, the volumetric flow rate from vessels increases, facilitating radiopharmaceutical transport to the tumor tissue. According to Fig. 12c,d, although it appears that the mean intracellular uptake of FDG increased with decreasing tumor size, this is not always fixed for all cases in the clinical conditions and requires consideration of direct effect of cell microscopic variables, such as GLUT expression and HXK enzymes. Consideration of these variables can be used as an advanced model in future efforts to more highlight FDG uptake. The proposed comprehensive spatiotemporal distribution framework enables the investigation of the relationship between molecular processes and imaging data. Inverse problems can be subsequently developed (future work) to estimate the parameters of interest, including kinetic parameters (K i ) from dynamic concentration data, diffusion coefficient in the tissue and the input function, towards personalized assessment of disease. Different imaging scenerios (tasks, exam durations, scanners, noise levels, etc.) can be simulated. Subsequently, by constructing appropriate numerical inverse methods (e.g., non-linear regression paradigms), the parameters of interest can be determined. This is hypothesized to result in improved estimates given enhanced SDM modeling over conventional kinetic compartmental methods. Future efforts also include artificial intelligence (AI) towards faster computation of the forward model, and improved inverse models. Such AI-assisted models can be developed via access to large real and synthetic PET datasets.
Overall, combining both the paradigms of therapeutics and diagnostics (i.e., theranostics) into clinical formulations allows for rapid assessment and adjustment of treatment to individual needs for personalized medicine. In theranostics, a diagnostic radiopharmaceutical (e.g. 68 Ga-PSMA-11 for PET imaging) is used to detect a www.nature.com/scientificreports/ molecular target with sufficient uptake and abnormal specificity. A therapeutic radiopharmaceutical (e.g., 177 Lu-PSMA-617) is then administered at a therapeutic dose level to treat the abnormal tissue. Quantitative imaging of the therapeutic radionuclide to image pharmacodynamic response to treatment allows the patient to select personalized and patient-specific dosages to reduce radiation toxicity and optimize the benefits of radiation therapy. Existing radiopharmaceutical therapy practice tends to focus on fixed doses for patients (e.g. 7.4 GBq (200 mCi) for 177 Lu-DOTATATE). However, to truly personalize therapies it will be necessary to predict the delivered dose. This will require construction of sophisticated digital twins for patients, incorporating more detailed biology and pharmacokinetic modeling. As such, the construction of models such as ours can help enable such a vision, such as our recent efforts [83][84][85] . While our current model still requires a significant amount of patient-specific data for clinical utility, this study shed lights on the field in what we believe is a needed angle toward a future where deeper understanding on effective factors and parameters translates to more realistic models and decisions. Thus, considering the outcomes of such models, more thorough models incorporating individual patient data will be provided.
Validation of numerical model. The distribution trend of FDG uptake in the tissue for C e and C m agrees well with the experimental results of studies by Carson et al. 86 and Eastman et al. 87 . Figure 13 shows the comparison of the mean total concentration of FDG in the tumor area with the experimental observation of Backes et al. 81 . Total FDG concentration is calculated in both extracellular and intracellular spaces, as measured in radionuclide imaging. There is a good agreement between the present study and the experimental observations, as the trend of both curves is nearly identical. However, since the domain and conditions of the experimental and modeling are different, the results are not exactly matched. After 15 min, the total concentration of FDG for the current study is very close to the total concentration in the experimental results, and after 25 min, the results exactly matched the experimental results. It should be noted that accurate examination of FDG concentration peaks requires accurate values of the input function at all times. Since the input function of Backes et al. is only available for up to one hour, it is possible to investigate the model's capabilities in distinguishing the kinetic steps of [ 18 F]FDG delivery to the extracellular matrix, its transport from the extracellular to the intracellular space, www.nature.com/scientificreports/ and its intracellular phosphorylation. Additionally, it is possible to evaluate spatiotemporal distribution of radiopharmaceutical by considering different heterogeneity of tumor microenvironment.  www.nature.com/scientificreports/

Conclusion
We have presented a comprehensive mathematical model to accurately simulate the spatiotemporal distribution of the [ 18 F]FDG radiopharmaceutical uptake within vascularized solid tumors. This model includes the production of tumor-induced microvascular network by a two-dimensional hybrid sprouting angiogenesis method along with the calculation of blood flow in the adaptable new vessels that remodels the network structure in response to metabolic and hemodynamic stimuli. Biological validation of the produced angiogenic networks was performed with in vivo observations. The spatiotemporal distribution of FDG was obtained by coupling intravascular blood flow and interstitial fluid flow equations with the SDM solute transport mathematical model. This model includes convection and diffusion transport mechanisms from the vessels to the tissue and within the tissue, as well as the reaction mechanism. We demonstrated the capabilities of the model by examining the effects of tumor size and MVD on the radiopharmaceutical distribution over time. Results show that the behavior of FDG radiopharmaceutical distribution is highly dependent on MVD. Due to the high permeability of tumor capillaries, the maximum total concentration occurred in the tumor area at all time frames. Additionally, the maximum concentration value occurs in the vicinity of the vessel walls and decreases with increasing distance from the capillaries due to the low diffusion coefficient. As the MVD increases, the mean FDG uptake into the tumor is also enhanced. In addition, this model correctly predicts the effect of desmoplasia, which involves an increase in IFP in larger tumor sizes. At smaller tumor sizes, an increase in volumetric flow rate, due to lower IFP values, and an increase in the contribution of diffusion mechanism from vessels increased the mean total FDG concentration. The proposed spatiotemporal modeling framework can be used to comprehensively evaluate the effect of various parameters on the spatiotemporal distribution of different radiopharmaceuticals, beyond the usage of conventional methods, including ODE-based kinetic compartment modeling. A goal of the present study is to help bridge the gap among different related fields, including mathematical oncology, radiology, and cancer diagnosis.