Investigation of the Spatiotemporal Responses of Nanoparticles in Tumor Tissues with a Small-Scale Mathematical Model

The transport and accumulation of anticancer nanodrugs in tumor tissues are affected by many factors including particle properties, vascular density and leakiness, and interstitial diffusivity. It is important to understand the effects of these factors on the detailed drug distribution in the entire tumor for an effective treatment. In this study, we developed a small-scale mathematical model to systematically study the spatiotemporal responses and accumulative exposures of macromolecular carriers in localized tumor tissues. We chose various dextrans as model carriers and studied the effects of vascular density, permeability, diffusivity, and half-life of dextrans on their spatiotemporal concentration responses and accumulative exposure distribution to tumor cells. The relevant biological parameters were obtained from experimental results previously reported by the Dreher group. The area under concentration-time response curve (AUC) quantified the extent of tissue exposure to a drug and therefore was considered more reliable in assessing the extent of the overall drug exposure than individual concentrations. The results showed that 1) a small macromolecule can penetrate deep into the tumor interstitium and produce a uniform but low spatial distribution of AUC; 2) large macromolecules produce high AUC in the perivascular region, but low AUC in the distal region away from vessels; 3) medium-sized macromolecules produce a relatively uniform and high AUC in the tumor interstitium between two vessels; 4) enhancement of permeability can elevate the level of AUC, but have little effect on its uniformity while enhancement of diffusivity is able to raise the level of AUC and improve its uniformity; 5) a longer half-life can produce a deeper penetration and a higher level of AUC distribution. The numerical results indicate that a long half-life carrier in plasma and a high interstitial diffusivity are the key factors to produce a high and relatively uniform spatial AUC distribution in the interstitium.


Supporting Information Numerical methods
To compute Equation (6), under the geometric configuration shown in Figure 1, in a conventional finitedifference or finite-element way requests a structured or unstructured body-fitted mesh generation at first. This pre-processing is not a trivial job and usually takes a serious amount of man and computer hours to do so especially when a good-quality mesh is demanded. After discretization, Equation (6) would end up with a large sparse linear system, which can be solved directly by some sparse solvers. However, this approach would request huge computer memory. An alternative is to solve the linear system by iteration methods, which do not demand large computer memory, but sometimes suffer poor convergence due to ill-conditioning. The convergence usually can be improved by pre-conditioning, but it takes skills to choose a good pre-conditioner and the associated iteration method. To avoid all the difficulties above, here we embed Equation (6) into a time-dependent problem as follows, and compute Equation (S1) until the steady-state is reached. Method of lines (MOL) is employed to compute Equation (S1). The basic idea of MOL is first to semi-discretize a time-dependent problem in space. This will end up with a set of coupled ordinary differential equations (ODE), which can be efficiently solved by many well-developed ODE solvers.
Here we applied the highly accurate multi-block Chebyshev pseudospectral method to discretize Equation (S1) in space. We first divided the physical domain shown in Figure 1 into three parts, denoted by sub-domains I, II, and III, as shown in Figure 2A. Since Chebyshev pseudospectral method can only be implemented on a square domain [−1, 1] × [−1, 1], we need to employ coordinate transformation to each physical sub-domain, denoted by (x j , y j ), j ∈ {I, II, III}, to rectangular computational domains, denoted by (ξ j , η j ), with −1 ≤ ξ j , η j ≤ 1, j ∈ {I, II, III}, respectively. Sigma transformation, frequently used in hydraulics, is applied to transform sub-domains I and III to rectangular ones.
Chebyshev pseudospectral method with Gauss-Lobatto grids (CPGL) has been a popular high-order numerical method for spatial discretization of partial differential equations [1]. To illustrate the simple idea of CPGL together with MOL, here we will use one-dimensional heat conduction problem as an example Note that when Equation (7) or Equation (S1) is considered, w will correspond to C or P i . In basic theory of CPGL, if w N (x) is the projection of w(x) into a polynomial space with degree less than or equal to N through Chebyshev polynomials as basis functions, then w(x) can be well approximated by where L N,k (x) is the k-th Lagrange interpolation polynomial of degree N , with the interpolation node x k being Chebyshev Gauss-Lobatto collocation points, Similarly, w ′ N (x) and w ′′ N (x) will be good approximations of w ′ (x) and w ′′ (x), respectively, and can be obtained through differentiation of Equation (S3), and where ′ and ′′ denote the first and second derivatives. In the spirit of collocation method, we are only interested in function values on those collocation points. Therefore, and where L ′ N,k (x l ) and L ′′ N,k (x l ) are so-called 1st and 2nd order collocation derivative matrices usually denoted as D (1) and D (2) . The entries of D (1) and D (2) can be found in Ref. [2]. In fact, D (2) can be well approximated by D (1) D (1) . Equations (7) and (S1) will further turn into a set of differential-algebraic equations (DAE) with the algebraic equation part from boundary conditions which are independent of time. This DAE system can be solved by many well-developed DAE solvers, and here we employed ode15s, a popular MATLAB script, to solve this DAE system [3]. Besides boundary conditions mentioned above, extra interface conditions between sub-domains rising from domain decomposition are also required. They are simply continuity of P i , C, and their normal derivatives. The current MOL+CPGL approach did show a great efficiency in computing Equations (7) and (S1), since it does not request huge number of grids to gain enough accuracy. Also Gauss-Lobatto collocation mesh tends to cluster near boundary and interface, as shown in Figure 2A. This may also help to resolve possible boundary layer economically.

Auxiliary hydrodynamic model
In order to validate our numerical method, we further developed the hydrodynamic model in the interstitial region. From Starling law, the pressure drops across arterial and venous walls are proportional to the transmural volume flow rate as and where subscripts "a" and "v" are employed to indicate the capillaries at arterial and venous ends, respectively. The capillary osmotic pressures are assumed to be π a = π v ≃ iCRT , with i and C denoting the van't Hoff factor and drug carrier concentration in capillaries, respectively. While the capillary pressure P a or P v can be measured in practice, the interstitial hydrostatic pressures P ′ a or P ′ v in the tumor interstitium right adjacent to vessel walls of arterial and venous microvessels respectively are yet to be determined.
To determine the values of P ′ a − P ′ v , we carried out intensive simulations to obtain the interstitial flow rate (Q tissue ) as an empirical function of the interstitial pressure drop (P ′ a − P ′ v ), hydraulic conductivity (K), capillary radius (r), and vascular distance (G), which was found to be with R 2 = 0.9639. Given the vascular distribution in our simulation model, the volume flow rates across the arterial capillary wall, interstitial space and venous capillary wall are equal to each other, i.e., By summing Equations (S8)-(S10), it yields the pressure difference between arterial and venous capillaries yields (S12) Alternatively, the transmural and interstitial flow rate is described by Once Q is determined by Equation (S13), P ′ a and P ′ v can then be determined from Equations (S8) and (S9) respectively, which concludes the Dirichlet boundary conditions setup for P i when solving Equation (S6) or Equation (S1) alternatively. From knowledge of normal tissue parameters listed in Table S1, we validated Equation (S13) by the computed volumetric flow rate Q = 8.08 × 10 −7 µm 3 /s, which is comparable to the value of 9.8 × 10 −7 µm 3 /s reported in Ref. [4]. This demonstrates the validity of our mathematical model and numerical method. Consequently, our method can help to identify the contributions of some key factors in the drug delivery and may facilitate the design of therapy strategies for tumors. Table   Table S1. The physiological parameters of the normal tissue used to validate our numerical method.