Simulation of Multispecies Desmoplastic Cancer Growth via a Fully Adaptive Non-linear Full Multigrid Algorithm

A fully adaptive non-linear full multigrid (FMG) algorithm is implemented to computationally simulate a model of multispecies desmoplastic tumor growth in three spatial dimensions. The algorithm solves a thermodynamic mixture model employing a diffuse interface approach with Cahn-Hilliard-type fourth-order equations that are coupled, non-linear, and numerically stiff. The tumor model includes extracellular matrix (ECM) as a major component with elastic energy contribution in its chemical potential term. Blood and lymphatic vasculatures are simulated via continuum representations. The model employs advection-reaction-diffusion partial differential equations (PDEs) for the cell, ECM, and vascular components, and reaction-diffusion PDEs for the elements diffusing from the vessels. This study provides the details of the numerical solution obtained by applying the fully adaptive non-linear FMG algorithm with finite difference method to solve this complex system of PDEs. The results indicate that this type of computational model can simulate the extracellular matrix-rich desmoplastic tumor microenvironment typical of fibrotic tumors, such as pancreatic adenocarcinoma.


INTRODUCTION
The process of cancer progression is driven by the communication between tumor cells and their surroundings. A dynamic tumor microenvironment typically consists of highly proliferating neoplastic cells of different phenotypes, necrotic tumor cells, infiltrating innate and adaptive immune inflammatory cells, cancer-associated fibroblasts, cancer stem cells, extracellular matrix (ECM), blood and lymphatic vessels, pericytes, healthy host cells, and a variety of soluble molecules (Hanahan and Weinberg, 2000;de Visser and Coussens, 2006;Tlsty and Coussens, 2006;Whiteside, 2008;Perez-Moreno, 2009). These cellular and molecular elements dictate the tumor progress from unregulated neoplastic growth to potential metastasis. In its heterogeneous milieu with complex tumor-induced interactions and mechanical stress, an evolving tumor mass also undergoes transient morphological changes arising from cell motility and cell-cell/cell-ECM interactions. Mathematical modeling of cancer progression including its associated microenvironment may be a useful tool for predicting tumor dynamics and cancer response to therapy.
To study the desmoplastic tumor microenvironment, we recently presented a tumor model (Ng and Frieboes, 2017) as a continuum scale multicomponent-multispecies system consisting of heterogeneous cell types and ECM. This thermodynamic mixture model, inspired by the one derived in Wise et al. (2008) and Frieboes et al. (2010), includes metabolic reactions, tumorigenic factors, desmoplastic response, as well as tumorinduced angiogenesis and lymphangiogenesis. The diffuse interface approach is implemented as derived in Wise et al. (2008), where thermodynamically consistent Darcy velocities and Fickian diffusive terms are determined from the energy variation. In the Helmholtz free energy equation, the square gradient model (Cahn and Hilliard, 1958;Cahn, 1959;Yang et al., 1976;Rowlinson, 1979) is used to describe interfaces arising from the adhesive properties of cells and ECM components, and an elastic energy term by Leo et al. (1998) is added to represent the elastic properties of the ECM component. Continuous blood and lymphatic vessel densities are modified from cell fluxes employed in Anderson and Chaplain (1998), Chaplain (1996), and Mantzaris et al. (2004). Sprout initiation conditions of vessels, as well as interactions between angiogenic factors, proteolytic enzymes, and the ECM component are inspired by Levine et al. (2000), Levine et al. (2001a,b). Lymphangiogenesis is assumed to behave similarly to angiogenesis. Nutrients and waste products from cell metabolism are governed by fluxes and reaction rates modified from Casciari et al. (1992).
This paper presents the fully adaptive non-linear FMG algorithm and numerical solution of the model presented in Ng and Frieboes (2017). This diffuse interface model is characterized by non-linear fourth-order Cahn-Hilliard type PDEs and narrow transition layers. A fully adaptive block-structure Cartesian mesh is used on adaptively refined grid levels to address the need of transient locally-refined mesh regions for the narrow transition layers. To avoid severe time-step restrictions of explicit methods and fine resolution of transition layers, a fully adaptive non-linear finite difference multigrid method inspired by the pioneering work of Wise et al. (2007) and Wise et al. (2011) is implemented to solve the system of equations. The equations are discretized in time by the semi-implicit Crank-Nicolson method. The error smoothing steps employ the non-linear Full Approximation Scheme, and iterations are carried out in V-cycles. Modified from the numerical solution and method given in Wise et al. (2007Wise et al. ( , 2011 we implement a fully adaptive non-linear full multigrid (FMG) algorithm with finite difference method to solve this complex set of PDEs.

Desmoplastic Tumor Model
Consider a tumor growing in a tissue domain ⊂ R 3 where tumor and healthy cells, as well as the ECM, are tracked with continuous volume fractions. In the diffuse interface model, adhesive forces hold the tumor cells together, creating a boundary layer of finite thickness between the tumoral and healthy regions. An appropriate distribution of ECM across the domain can be achieved by a carefully chosen free energy term. Dimensionless governing equations for multispecies tumor growth systems are presented here.
We assume that the liquid (extracellular fluid) volume fraction stays constant, as well as the total solid (ECM, tumor and healthy cells) volume fraction, and set the densities of all components to unity. Dimensionless variablesφ V ,φ D ,φ E , and φ H are the normalized volume fraction (normalized by the total volume fraction) of viable tumor cells, dead tumor cells, ECM, and healthy host cells, respectively. Their transient diffusionconvection governing equations are given as Ng and Frieboes (2017):

Cells and ECM Components
The total tumor volume fraction isφ T =φ V +φ D andM i =Mφ i is the positive non-constant mobility of component i. Chemical potential terms are given below: whereẼ is the infinitesimal strain,ũ d is the displacement vector, E * E andẼ * C are the Eigenstrain tensors for ECM and cells, respectively. In Equations (2.10) and (2.14), δ mn = 1 for m = n and δ mn = 0 for m = n.L E 1 ,L E 2 ,L C 1 , andL C 2 are Lamé constants for the ECM and cell components. The displacement vectorũ d is solved by setting ∇ ·T i = 0, whereT i = T 1iT2iT3i and i = 1, 2, 3.
The dimensionless cell-ECM phase pressurep and interstitial fluid phase pressureq, with their corresponding velocitiesũ α and u β are computed using the following equations:

Tumorigenic Species
Here,F E reflects the concentration within the ECM phase, assuming the volume of myofibroblastic cells is negligible and the ECM phase is continuous throughout the domain. Similarly, dimensionless ECM based concentrations of blood and lymphatic vessels,B E n andL E n respectively, with their corresponding diffusive flux terms are given by. Tissue effective diffusivitiesD n ,D g ,D w ,D ,D b ,D a ,D s ,D tgf , D taf ,D m ,D FE ,D BnE ,D LnE , and tissue effective mass transfer coefficientsλ B,n ,λ B,g ,λ B,w ,λ B, , are represented byψ and determined as follow:

Blood and Lymphatic Vessels
whereψ E ,ψ T , andψ H are dimensionless diffusivity or transfer coefficient in the ECM, tumor, and healthy-host cell domain respectively. The total cell volume fraction isφ C =φ T + φ H . The full list of parameters, constants, factors used in nondimensionalization, source terms, rate terms, and adjustment factors can be found in Supplementary Tables 1-6. We define the following Neumann boundary conditions for the cell and ECM volume fractions, and Dirichlet boundary conditions for solid cell pressure and chemical potentials at all external boundaries: where n is the outward normal of a boundary. For nutrients and waste products, as well as tumorigenic species, Dirichlet boundary conditions are imposed, with the exception of myofibroblastic cell species, where Neumann boundary conditions are applied: Blood and lymphatic vessels are assumed to be at their corresponding far-field values at external boundaries: whereB ∞ =L ∞ = 0.2 is used here.

Numerical Methods
The mathematical model is first discretized in time and space. After reorganizing the discretized equations, a non-linear relaxation procedure is employed.

Computational Domain
We consider a rectangular 3D domain = (L x , R x ) × L y , R y × (L z , R z ). Let the domain be discretized into N x × N y × N z cells. Bounded by boundary cells, the domain is covered by the following sets of cell centers, cell edge points, and cell grid , (2.2.6) and if partitions in the three directions are uniform and equal: for 1 ≤ i ≤ N x , 1 ≤ j ≤ N y , and 1 ≤ k ≤ N z , we get the following coordinates for cell centers and cell edges or corners defined in Equations (2.2.1)-(2.2.5): following their corresponding ranges indicated in Equations (2.2.1)-(2.2.5). The differential, Laplacian and flux terms associated with the model are in Supplementary Materials.

Model Discretization
The model consists of a set of stiff differential equations that are fourth-order in space. At time step a with time step size θ , they are discretized in time using the Crank-Nicolson Method as in Wise et al. (2011). Details of the discretization as well as the multigrid V-cycle iterations and non-linear Gauss-Seidel relaxations are given in Supplementary Materials.
The tumor model solved using the adaptive full multigrid V-cycle was coded in C and simulations were performed on a FIGURE 1 | For simplification and ease of visualization, a two dimensional example is given here. The multigrid algorithm consists of a hierarchy of grid levels κ = −2, . . . , 2, shown here on the top rowwith their corresponding subdomains/blocks. Level κ = 0 is the root level and levels κ = 1 and 2 are refinement levels. Global levels κ = 0, −1, and −2 cover the entire computational domain, whereas refinement levels κ = 1 and 2 cover subdomains where refinement is needed. Shown in the example here, the locally refined block-structured κ = 1 consists of one block, B 1,1 , and κ = 2 consists of two blocks, B 2,1 and B 2,2 . The middle row shows the corresponding meshes for each level. For refinement levels in the example here, mesh G 1,1 for block B 1,1 , G 2,1 and G 2,2 for blocks B 2,1 and B 2,2 . The composite grids of all levels is shown in the last row. node equipped with 768 GB of RAM and 32 Intel Xeon 3.3 GHz cores running CentOS 6.7 x86_64. The algorithms were partially parallelized using OpenMP to achieve higher performance.

Self-Adaptive Full Multigrid Full Approximation Scheme V-Cycle Algorithm
The self-adaptive full multigrid process combines FMG with the self-adaptive approach. The process involves constructing a multilevel block-structured mesh and seeking the solution on the new mesh structure. Mesh construction and refinement, as well as FAS multigrid in an adaptive FMG setting were developed as follows.

Multilevel Mesh Refinement
We start with a rectangular computational domain as described in Methods. Each level κ covers a domain of κ with mesh size η κ . Grid levels are numbered as κ = κ min , . . . , 0, . . . , κ max , where κ min represents the coarsest mesh level and κ max the finest. Global grid levels, κ = 0 − κ min , cover the entire computational domain κ = .
After solution is obtained for the current finest level, we use the undivided gradient test  to tag cells for refinement. The undivided gradient test is used to capture the diffuse interface region (tumor-host cell domain). Let the critical value of some indicative variable ψ κ x, y, z be C κ , cell-centered coordinates on level κ be x i,j,k , the following test is performed on cells on level κ: where F κ are coordinates of flagged cells. Another criteria that can be used to tag cells for refinement is the relative truncation error test (Trottenberg et al., 2001). In the relative truncation error test, the relative truncation error with respect to κ and κ−1 , as used in Equation (3.1.20) in the FAS cycle, is used to flag cells for refinement. A subroutine is created and called to perform this task: where the current finest-grid data is used and F a,n κ is a list of flagged coordinates.
The flagged cells are then rearranged into patches of rectangular refined mesh in the next finer level using a clustering algorithm from Bell et al. (1994) and Berger and Rigoutsos (1991) with minor modification. Create a list of blocks and start with just one block containing all flagged cells given by F a,n κ . The following procedure is used to divide flagged cells into blocks: (1) Compute efficiency of the block, which is the ratio of flagged cells to total cells. If the efficiency is below the threshold efficiency, and, if the size of the block is more than twice the threshold size, then continue. Else, accept and add the current block to the list, return to Step (1) for the next block on the list. (2) Compute signatures, which is the number of flagged cells in each slice along each direction. Checking for gaps along all dimensions starting from the longest edge, find a gap closest to the center of that direction. If gaps exist and an optimum gap location is found, slice the block into two along the gap, and go to Step (5). If no gaps are found, continue.
(3) Compute second derivatives of the signatures along all directions. Checking for inflection points starting from the longest edge, find an inflection point closest to the center of that direction. If inflection points are found and an optimum inflection location is found, slice the block into two along the inflection point, and go to Step (5). If no inflection points are found, continue. (4) If the efficiency of the block is above a minimum allowed, accept and add the current block to the list, returning to Step (1) for the next block on the list. Else, divide the block into two along the mid-point of the longest dimension and continue. (5) Delete all empty slices along the edges of each block, and add the two trimmed blocks to the list. Repeat from Step (1) until all blocks created are checked.
The above task is assigned to the following subroutine of block generation: where B κ+1 is an array of blocks on the refined level κ + 1, with rows of coordinates corresponding to corners of each block.
Following the generation of new refined grids on level κ + 1, the new grids are populated with data from level κ and the old level κ+1. All cell-centered data for variables on the newly refined grids must be generated by higher order interpolation, such as cubic interpolation given by Equations (3.1.28) and (3.1.29), from the coarse grids below. If a new level κ + 1 grid cell overlaps any old level κ + 1 grid cell, the previous time-step data for the new grid cell is copied from the old. However, for any new level κ + 1 grid cell that does not overlap with an old κ + 1 grid cell, its previous time-step data has to be obtained from the coarse grid below. For all old level κ + 1 grid cells that do not overlap with any new κ + 1 grid cells, data are averaged and stored in the coarse level κ grid. Refer to Figure 2 for illustration.
After populating the newly refined grid, ghost cells surrounding each grid patch must be constructed. We first compute data for the ghost cells using the quadratic interpolation given by Colella et al. (2009). Illustrated in Figure 3, quadratic interpolation is performed twice. First, to get data from the coarse-grid for the intermediate points a and b: respectively. Then, to compute fine-grid boundary data on ghost cell c, perform quadratic interpolation using points a, e, and g. Similarly, interpolate using points b, f , and h to obtain data at ghost cell d. Following Wise et al. (2007), all ghost cells of each fine-grid patch are first filled using the quadratic interpolation above. If a ghost cell falls on a cell on any neighboring fine-grid patches of the same level, then the ghost cell data is replaced by the more accurate cell-centered data from the cell of the neighboring patch.

Adaptive FAS V-Cycle
Consider a case with κ grid levels with κ = 0 being the root level. The levels κ = κ max and κ = κ min represent the level with the finest and coarsest grid, respectively. Each level of κ > 0 contains one or more blocks of refined rectangular grids covering a subdomain κ , whereas levels κ ≤ 0 cover the entire computational domain κ = .
Averaging is used for the restriction operator R κ−1 κ : which reduces to a four-point average in a 2D case. Linear interpolation is used for the prolongation operators P κ κ−1 in error correction steps within a V-cycle. In the 2D cellcentered discretization cases as depicted in Figure 4, the bilinear interpolation operators are given by Trottenberg et al. (2001): for points marked a, b, c, and d, respectively. For 3D cases as shown in Figure 5, trilinear interpolation produces the following FIGURE 2 | Populating a newly refined grid. For simplification and ease of visualization, a two dimensional example is illustrated here. An existing block, with its boundary marked by --, is overlapped by a newly refined block with boundary marked by . All cell-centered variable data on the newly refined grid are obtained via cubic interpolation from the coarse grid data below. The previous time-step data for the new grid points (•), which overlap old grid points, are copied directly from the old grid points they overlap; for new grid points ( ) that do not overlap with any old grid, their previous time-step data are obtained from the coarse grid below. The remaining data on the old grid (•) that do not overlap any of the new grid points are averaged and copied to the coarse grid data points (×). Frontiers in Physiology | www.frontiersin.org for points a and c, b and d, e and g, f and h, respectively. The coarse grid operator L κ−1 is generated in a way that is analogous to the computation of L κ on the fine grid. Alternatively, the Galerkin coarse grid operator can be used: where R κ−1 κ and P κ κ−1 are appropriately chosen transfer operators (Trottenberg et al., 2001).
A FAS V-cycle consists of two iterating components. The outer time iteration travels through one V-loop, starting from the finest mesh level, looking for the fixed point solution of L ψ a,r i,j,k , ψ a,r−1 i,j,k = R ψ a−1 i,j,k , ψ a,r−1 i,j,k as derived in Supplement (Multigrid V-Cycle Iterations). Starting from ψ a,0 , each successive V-cycle produces an approximation that is converging toward the fixed point solution, ψ a , at the current time step a. After r * outer time iterations (after r * V-loops), if the approximation ψ a,r * results in an error within a tolerable bound, an approximated solution for the current time step is established by ψ a = ψ a,r * . The inner iteration is the nonlinear Gauss-Seidel relaxation as described in Supplement (Nonlinear Gauss-Seidel Relaxations) and given in Supplementary Equation (1.4.39). This procedure occurs twice on each level. Relaxation first takes place during the error smoothing step while traveling down from the finest to the coarsest grid level, and again after coarse grid correction while going up from the coarsest to the finest grid level. Hence, the starting approximation for the error smoothing step is obtained from the previous iteration ψ a,r,0 κ =ψ a,r−1 κ or restricted from the relaxed fine-grid approximation ψ a,r,0 κ =R κ−1 κ ψ a,r,ν 1 κ+1 , whereas the coarsegrid-corrected approximation ψ a,r,0 κ =ψ a,r,CGC κ is used during the correction step.
The computation of a new time iterate ψ a,r κ using an adaptive FAS multigrid-cycle, starting with any of the refined or global grid levels described in the Multilevel Mesh Refinement subsection, can be done recursively using a two-level model as summarized in the operator below (Brandt, 1977;Trottenberg et al., 2001;Wise et al., 2007).

III. Post-smoothing:
Compute ψ a,r κ by applying ν 2 smoothing steps to ψ a,r,CGC κ on κ : -Update the ghost cells on κ level using interpolation and exchange for ψ a,r κ on neighboring patches. V-or W-cycle is defined via the cycle index γ = 1 and γ = 2, respectively. The initial value used in V-cycles may come from the converged solution at previous time ψ a,0 = ψ a−1 or from coarse-grid approximations as in FMG method described in the next section. The SMOOTH routine in Equations (3.1.18), (3.1.22), and (3.1.27) uses Red-Black (or odd-even) ordering for the relaxation process, with ν b extra smoothing steps for d κ cells within the boundaries.

Adaptive FMG Structure
In FMG, a better initial approximation for finer-grid iteration is obtained from a coarser-grid level (Kronsjö and Dahlquist, 1972;Brandt, 1977). The procedure begins with a set of solutions on the coarsest grid, followed by interpolation of the solution set to a fine-grid level, providing an initial guess for the fine-grid multigrid approximation. After a few multigrid cycles on the finegrid level, the solution is again interpolated to the next fine-grid level. This proceeds until the finest grid level has been reached, if refinement is needed for each level κ > 0. Discretization accuracy can be reached on each grid level within a few multigrid cycles in this manner. The structure of full multigrid is illustrated in Figure 6.
A higher order scheme is normally employed for an FMG interpolation. Cubic Lagrangian interpolation is one such method used in the FMG interpolation steps to compute an approximation at the fine grid. In 1D as depicted in Figure 7, approximation ψ κ (x) at the fine grid point can be computed as for point a and b respectively. Cubic interpolation in the y-and z-directions for bicubic and tricubic interpolation, in 2D and 3D cases, respectively, are analogous to the 1D calculation shown. FMG can be performed on locally refined grids which are defined a priori. Here, we use self-adaptive FMG, where refined grids may or may not be generated, one layer at a time, based on some predefined criteria. Self-adaptive FMG advances as follows: (1) Perform FMG on a series of global grids until a satisfactory accuracy has been obtained on grid level κ = 0. (2) Determine if refining the finest grid would be beneficial, based on the criteria outlined in the Multilevel Mesh Refinement subsection. If no, iteration stops. If yes, refine grid in local regions. (3) Use cubic interpolation on the current fine-grid (κ) approximation to compute an initial guess for the refined grid level (κ + 1), proceed by performing a number of multigrid cycles on the refined level until the required accuracy has been satisfied. (4) If the finest level κ max has been reached, the time iteration stops. Otherwise, return to Step (2).

End If End For
A parameter σ is used to adapt the accuracy to the different level of meshes. Tolerance is reduced by σ for each finer mesh level, tol κ = tol κ−1 σ . The parameter σ and the initial value of tolerance have to be selected with care in order to ensure the overall accuracy of the final grid structure.

System Solution
Next, we applied the fully adaptive non-linear FMG algorithm to solve the desmoplastic tumor model. The set of tumor model related parameters, constants, source terms, rate expressions, and adjustment factors are listed in Supplementary Tables 1-6.
We assume that the tissue consists of a fixed fraction of liquid interstitial-fluid phase (φ β = 0.2) and solid cell-ECM phase (φ α = 0.8). We also assume that the solid cell-ECM phase is only composed of viable tumor cells (φ V ), dead tumor cells (φ D ), extracellular matrix ECM (φ E ), and healthy host cells (φ H ). Nutrients, waste products, and tumorigenic species are assumed to be carried in the liquid phase, while the myofibroblastic cell species (F E ), blood (B E n ) and lymphatic (L E n ) microvessels are assumed to take up negligible volume and exist within the continuous ECM component.
We assume that the source and sink of viable tumor cells comes from cell proliferation and losses to necrosis respectively. The dead tumor cell species hence comes from the necrosed viable tumor cells and undergoes lysis. The ECM is assumed to be secreted by the myofibroblastic cell species and degraded by the matrix degrading enzyme species (m). The healthy host cells species is assumed to be homeostatic. Nutrients such as oxygen (ñ) and glucose (g) enter the tissue via blood vessels and are consumed mainly by the viable tumor species. Tumor growth factors (tgf ) and tumor angiogenic factors (taf ) are assumed to be secreted by viable tumor cells and can undergo degradation. Tumor angiogenic factors also have an uptake term by proliferating vessels. Myofibroblastic cell species is assumed to go through proliferation, apoptosis, and necrosis, while the blood and lymphatic vessels undergo restructuring and potential loss due to crushing by the surrounding tissue pressure. To slightly perturb the symmetry of tumor progression, we let the fraction of blood vessels that are sprouting be different in different regions of the domain, ranging from 0.1 to 0.8.
As described in Methods, all PDEs were discretized using the Crank-Nicolson method. To solve the set of cell-centered discretized equations, we applied full multigrid algorithm and adaptive full approximation scheme V-cycle with Gauss-Seidel Red-Black smoothing to the multiple grid level system. Here, we use a total of five grid levels as depicted in Figure 1, with three increasingly refined global levels (κ = −2, −1, 0) covering the entire domain and two adaptively refined levels (κ = 1, 2). The domain is = (0, 40) 3 and the mesh sizes for the finest global level to the finest adaptively refined level are η 0 = 40/32, η 1 = 40/64, and η 2 = 40/128. The time step-size used is θ = 1 × 10 −2 and minimum tolerance is set at 5 × 10 −4 . Key solver related parameters are listed in Table 1.
The center of the computational domain is simply seeded with viable tumor cells (φ V = 0.65) at the beginning of the simulation. We assume that there are no dead tumor cells (φ D = 0) initially and the ECM is distributed evenly (φ E = 0.35) across the domain. The healthy host cells thus take up the remaining volume in the solid phase of the tissue. The initial tumor shape is therefore rectangular with a sharp interface. The simplified source term for the viable tumor cell species used in Equation (2.1), where only mitosis and necrosis are considered, is given by: where H is the Heaviside function,ñ h = 0.3 is the hypoxic level of oxygen,λ M,V = 1 andλ N,V = 3 are the rate constants of mitosis and necrosis, respectively, for viable tumor cells.ñ v,V = 0.21 andg v,V = 0.1 are the oxygen and glucose viability limits, respectively, for the viable tumor cell species. The mitosis rate is assumed to be upregulated by the level of tumor growth factors. The center of the tumor mass experiences a significant drop in nutrient levels. When the oxygen level drops below the hypoxic level, viable tumor cells cease to reproduce. Necrosis takes place if the oxygen level falls below the viability limit, where viable tumor cells necrosed to dead tumor cells. The source term for the dead tumor cell species in Equation (2.2) is hence given bỹ whereλ L,D = 1 is the lysis rate constant for the dead tumor cell species. The lysed dead tumor cells are assumed loss to the interstitial fluid. The other major solid component of the tissue is the ECM, and its source term used in Equation (2.3) is whereλ F,E = 5 is the secretion rate constant of ECM by the myofibroblastic cell speciesF andλ de,E = 5 is the degradation rate constant of ECM by the matrix degrading enzyme species m. Factors affecting the rate constant of ECM secretion by the myofibroblastic cell species are included in the adjustment factor A F,E given by where F F n,E = 2 is the effective factor of hypoxia on upregulating the production of ECM by myofibroblastic cells andñ v,F = 0.21 is the oxygen viability limit of the myofibroblastic cell species. The constanttgf F,E = 0.2 is a lower threshold of tumor growth factors, below which the production of ECM by the myofibroblastic cell species is assumed negligible. A complete list of source terms, rate expressions, and adjustment factors is shown in Supplementary Table 6.
The transient tumor progression is shown in Figure 8, via tracking the evolution of the tumor isosurfaceφ T = 0.15. With the specific set of constants and parameters used in this case, the tumor undergoes regression in the time frame shown. Nutrients levels drop below the viability limit in the center of the tumor mass, causing dead tumor cells to accumulate from necrosed viable tumor cells, which are then lysed from the tumor mass. The increase in tumor growth factors in the tumor region upregulates the secretion of ECM by the myofibroblastic cell species, resulting in a high ECM environment within the tumor mass (as in Figure 9).
As mentioned in Results (Adaptive FMG Structure), a higher order scheme is normally used in the FMG interpolation step when interpolating the solution to a finer grid level. Trilinear interpolation as shown in Equations (3.1.12-3.1.15) and cubic interpolation given in Equations (3.1.28) and (3.1.29) are used in the FMG interpolation step and the simulation results for time = 10 are compared in Figure 9. The ECM and viable tumor , forκ > 0 Number of cells inward from a boundary to be included in the near-boundary extra smoothing steps Set to 2 κ−κ min for κ ≤ 0 Set to 1 8 of the minimum block dimension (in x-, y-, and z-direction) for κ > 0 r κ = 1 Number of V-cycles (iterations) performed before the finest global level is reached (for levels κ min < κ < 0, which is r −1 here) σ = 4 Factor by which the tolerance is reduced from level κ to κ + 1 C κ = 0.05 Critical value used in the undivided gradient test in the FLAG routine threshold_eff = 0.9 Upper limit of efficiency, any blocks below this limit will go through the splitting algorithm. (Used in the BLOCKGEN routine) threshold_size = 10 3 Minimum size of a block allowed. (Used in the BLOCKGEN routine) min _eff = 0.5 Lowest limit of efficiency allowed, any blocks below this limit will be bisected. (Used in the BLOCKGEN routine) levels are slightly higher in the case where cubic interpolation is used, resulting in higher solid cell-ECM phase pressure. Since the myofibroblastic cell species and vessels are assumed to reside within the ECM component, the higher ECM volume fraction in the high order scheme case also results in higher myofibroblastic cell species (presented in green as myF) and vessel densities (presented in red as Bn for blood and blue as Ln for lymphatic) in tissue. Concentration profiles of interstitial fluid based species show no visible difference in the two cases. The adaptive block-structured mesh system for the simulation in Figures 8, 9 at time = 50 is shown in Figure 10. For global level κ = 0, there is one block (shown in black) covering the entire domain with mesh size η 0 = 40/32. In the first level of refinement at level κ = 1, there are four blocks (shown in green) with mesh size η 1 = 40/64. In the next level of refinement at level κ = 2 and within the coarser domain 1 , there are six blocks (shown in red) with mesh size η 2 = 40/128. As shown in Table 1, the critical value of C κ = 0.05 in Equation (3.1.1) is used in the undivided gradient test to flag cells for refinement. We also buffer each flag cell by flagging four cells surrounding it in each direction, creating a cube of 9 × 9 × 9 flagged cells, or less if it is near any external or ghost boundaries. The thresholdsize sets the minimum block size to 10 3 , resulting in overall bigger blocks generated by the BLOCKGEN routine in Equation (3.1.4).
In Figure 11, the transient degree of freedom (DOF) for the simulation in Figures 8, 9 is plotted. The degree of freedom is represented by the total cell-centered grid points in the 3D domain. The mesh with just the three global levels (κ = −2, −1, 0) has only 37,376 DOF. In the first iteration, both refinement levels (κ = 1, 2) begin to contain flagged cells and the mesh now has 224,640 DOF. At time = 50, there are 554,568 DOF, increased by a factor of 2 from the first iteration.
A 3D convergence test as described in Wise et al. (2007) was done by varying grid and time-step sizes. Three simulations were performed where grid spacings used for the root level κ = 0 are 16 3 , 32 3 , and 64 3 , respectively. We focused on the convergence of the tumor volume fraction,φ T , and in all three cases, data were interpolated to the cell centers of their global FIGURE 9 | 1-D Profiles of species with trilinear (left side) vs. cubic (right side) interpolation used in FMG-interpolation at time = 10 (sliced at j = k = 58). First row: tumor viable speciesφ V , dead speciesφ D , and ECM speciesφ E . The overall tumor pressure is labeled by p α . Second row: diffusible substances driving the tumor evolution, including oxygen (O 2 ), glucose (Glu), carbon dioxide (CO 2 ), bicarbonate (Bic), lactate (Lac), and hydrogen ions (H + ). Third row: Concentration of myfibroblasts (myF), tumor growth factors (TGF), and matrix degrading enzymes (MDE). Fourth row: corresponding density of blood vasculature (Bn), lymphatic vasculature (Ln), and tumor angiogenic factors (TAF). The myofibrobrastic cell species density in tissue (shown in the third row of plots as myF) is computed bỹ F =F EφE . Similarly, the blood and lymphatic vessel densities in tissue (shown in the last row of plots as Bn and Ln) are given byB n =B E nφE andL n =L E nφE respectively.
uniform grids corresponding to their finest mesh η κ max . With κ max = 2, the grid sizes η 2 and the corresponding time step sizes θ , which follow the linear refinement path θ = 0.032η, are outlined in Table 2. The errors were calculated by comparing cell center data to the averaged cell center data of the finer grid set, and the l 2 norms of the errors were used to obtain the convergence rates as shown in Wise et al. (2007). The errors and the rate of convergence are also listed in Table 2. The results show that the algorithm presented herein is first-order accurate. Even though the discretization is second order in both time and space, the first-order rate of convergence attained, as also reported by Wise et al. (2011), is expected since the surface adhesion terms in Supplementary Equation (1.2.6) were treated explicitly. Non-smooth functions used in adjustment factors and first order interpolation function used in prolongation are other potential contributing factors to the first-order convergence.

DISCUSSION
This paper illustrates the application of a fully adaptive, nonlinear full multigrid, finite-difference algorithm to solve a diffuse interface desmoplastic tumor system (Ng and Frieboes, 2017). The set of PDEs in the model is discretized in time using the Crank-Nicolson method. A Non-linear Full Approximation Scheme is used in the full multigrid V-cycle iterations, and Red-Black ordering is used in the Gauss-Seidel relaxation. A blockstructure multilevel Cartesian mesh isused consisting of three global levels with mesh sizes η −2 = 40/8, η −1 = 40/16, η 0 = 40/32, and two adaptively refined levels with mesh sizes η 1 = 40/64, η 2 = 40/128. A numerical simulation of desmoplastic tumor progression in 3D which shows that the algorithm is capable of simulating an ECM-rich tumor environment and the handling of morphological changes during tumor progression (Ng and Frieboes, 2017).  The diffuse interface model is chosen to bypass the need to track the transient position of sharp interfaces, as well as enforcing complicated boundary conditions across these interfaces, such as those required in a sharp interface model. The adaptive multigrid algorithm used can efficiently handle the narrow transition layers as well as larger morphological evolutions. As opposed to lexicographic ordering, the Red-Black sweep used here in relaxations/error smoothings allows for easy parallelization to minimize the computational expense.
In the future, 3D meshes with more refinement levels and finer grids will be evaluated. Different criteria, such as the relative truncation error test or a simple volume fraction test, will be used to flag cells for refinement, and the minimum block size allowed during block generation will also be reduced. The mathematical model may be augmented to include additional tumor cell species as well as immune cell species, together with their corresponding functions and effects. Tumorigenic species such as hormonal growth factors and chemoattractants could be added. An anticancer drug species may be added to study various tumor responses to different therapeutics. A more expansive role and interaction of the lymphatic system with the microenvironment (e.g., Scianna et al., 2013;Swartz, 2014) could be implemented. The current numerical model is coded in C and is partially parallelized. Full parallelization of the program will be explored to speed up and extend the capabilities of the code to handle models of larger scale.

AUTHOR CONTRIBUTIONS
HF: study conception and scientific oversight; CN and HF: model definition and design; CN: initial implementation; CN and HF: code debugging; CN and HF: manuscript writing and revision. Both authors reviewed and approved the final manuscript.

FUNDING
This work was supported in part by the National Institutes of Health/National Cancer Institute (grant numbers U54CA143907, R01CA180149, and R15CA203605), and by the Department of Bioengineering and the James Graham Brown Cancer Center at the University of Louisville.