A multiphase model for three-dimensional tumor growth

With the aging world population, a surge in cancer incidence is anticipated in the coming years, with a major societal and economic impact. With such a scenario, the development of novel therapeutic strategies is critical for improving the prognosis, outcome of intervention and quality of life, and for minimizing the economic impact. In this context, computational models for tumor growth and its response to different therapeutic regimens play a pivotal role. Over the last two decades, multiple models have been developed to tackle this problem. As discussed in the comprehensive works of Roose et al (2007), Lowengrub et al (2010) and Deisboeck et al (2011), three major classes of models have been proposed: discrete, continuum and hybrid models. Discrete models follow the fate of a single cell, or a small cohort of cells, over time. As such, they cannot capture tissue mechanics aspects, nor are the modeled subdomains representative of the whole tumor. However, they explain cell-to-cell cross signaling and cell response to therapeutic molecules (Perfahl et al 2011). On the other hand, continuum models describe cancerous tissues as domains composed of multiple homogeneous fluids and solid phases interacting one with the other. Differential equations describe the spatiotemporal evolution of the system, but no direct information is provided at the single cell level (Roose et al 2007). Finally, hybrid models incorporate different aspects of discrete and continuum models, depending on the problem of interest. For instance, they represent cells individually and extracellular water as a continuum (Chaplain 2000, Anderson 2005, Bearer et al 2009).

At the very early stages, solid tumors are composed of a few abnormal cells growing within an otherwise healthy tissue. The vasculature is generally absent, and the tumor cells (TCs) take all their nutrients by diffusion from the surrounding tissue. This is defined as the avascular phase for a solid tumor. As the mass of TCs increases, the extracellular matrix (ECM) undergoes extensive rearrangements with increased deposition of collagen fibers, making the resulting tissue thicker and more difficult to trespass (Jain 1999, Jain and Stylianopoulos 2010). Also, since the TCs divide much faster than normal cells, the growing tumor mass exerts mechanical stresses on the surrounding healthy tissue, leading to the localized constriction and, at times, collapse of blood and lymphatic vessels. At this point, the TCs are already numbering millions and the malignant tissue has reached a characteristic size of hundreds of microns. A necrotic zone appears deep inside, far from the pre-existing vasculature and the interstitial fluid pressure (IFP) builds up against the vascular hydrostatic pressure, mainly due to the compression of the healthy tissue, obstruction of the lymphatic vessels and hyper-permeability of the new blood vessels. Using proper biochemical stimuli, the TCs recruit new blood vessels (angiogenesis) to support a continuous transport of nutrients and oxygen. This is defined as the vascular phase of a solid tumor. Over time, these new blood vessels also become a preferential route for the malignant mass to shed into circulation millions of abnormal cells that, transported by the blood flow, would reach distant sites and eventually lead to the development of secondary tumors. This is the metastatic phase, typically occurring for a few solid tumors. This briefly describes multiple phases and stages that characterize the evolution of tumors which cannot be accurately captured in a single, comprehensive computational model. Here, the focus will be on tumor initiation, and on a novel continuum model for the evolution of avascular tumors.

Most continuum models for avascular tumors describe the malignant mass as a homogeneous, viscous fluid and employ reaction–diffusion–advection equations to predict the distribution and transport of nutrients and cells (Roose et al 2007). Cell diffusion, convection and chemotactic motion are included, and cell proliferation is governed by mass and momentum balance equations. The first model was by Casciari et al (1992). More advanced models included intracellular mechanical interactions (pressure, shear and adhesion) and interaction of cells with the interstitial fluid (IF) pervading the ECM. In these cases, momentum balance equations and constitutive relations are also required for describing the stress–strain response of each individual phase. One of the earlier models (Byrne and Chaplain 1996) treated the TC as a viscous liquid and introduced, quite artificially, a hydrostatic pressure within the tumor domain representing the IFP. More sophisticated models since then have treated the solid and fluid phases independently. For instance, Roose et al (2003) modeled the tissue matrix as a linear poroelastic solid, while the IF was prescribed to obey Darcy's law. Cell growth was incorporated into the stress–strain relationship, still imposing small displacements. See also Sarntinoranont et al (2003).

Byrne et al (2003) have proposed a new class of models derived in the multi-phase framework of mixture theory. Mixture theory consists of a macroscopic description (level of observation) of the system where conservation laws are introduced in analogy with the balance laws of single bodies. Additional terms are introduced to account for the interaction among phases. The disadvantage of this approach is that no connection is made with the microscopic reality. Interfacial properties are absent from both conservation laws and constitutive equations—a serious deficiency when applied to porous media (Gray and Miller 2005). Within this approach the cellular phase (for both tumor and healthy tissues) is modeled as a viscous fluid and the IF as inviscid. Although the mixture theory formalism is potent and flexible, major challenges lie in the treatment of the interfaces arising from the different phases. Traditionally, two classes have been proposed: the sharp interface method, considering the interface as a sharp discontinuity; and the diffuse interface method, considering the interface as a diffuse zone. The sharp interface approach—difficult to implement for interfaces separating pure media (IF) and mixtures (TC and healthy tissue)—has been followed by Preziosi and Tosin (2009) and Preziosi and Vitale (2011a, b). However, necrotic cells are not distinguished from live TCs: tumors are modeled as if necrotic cells are no longer part of the tumor. They are hinted at in the source/sink term but the related balance equations are missing. Their inclusion would require accounting for an additional interface between living and dead cells, which is not sharp in nature. On the other hand, the diffuse interface approach introduces an artificial mixture at the interface, and the challenge here is to derive physically, mathematically and numerically consistent thermodynamic laws for these interfaces. Wise et al (2008), Cristini et al (2009), Oden et al (2010) and Hawkins-Daarud et al (2012) have all followed this approach. However, they include only one interface, separating the TC from the healthy tissue. Strictly, this is insufficient in the mixture theory formalism where each interface should be accounted for throughout the whole computational domain. The models lack some rigor because the designation of phases as distinct from chemical constituents comprising a phase is unclear. Consequently, some of the balance equations contain terms that cannot be justified on a theoretical basis. These simplified approaches lead to fourth-order-in-space parabolic partial differential equations of the Cahn-Hilliard type. This entails some difficulties for three-dimensional (3D) solutions with finite element methods because higher order basis functions are needed than in the realm of second-order spatial operators (Gomez et al 2008). Further, considering more than two or three phases becomes cumbersome, especially if a solid phase is included.

There is a need for tumor growth models for the dynamics of multiple phases and interfaces in a physically and numerically sound way. Recently, the thermodynamically constrained averaging theory (TCAT) framework was established (Gray and Miller 2005, Gray et al 2012) for continuum, porous media models that are thermodynamically consistent across scales. Here, the TCAT formalism will be used to predict the growth of tumors under different physiologically relevant conditions. We show that second-order differential equations can accommodate more phases than most of the existing models. The interface behavior is modeled through surface tension (Dunlop et al 2011, Ambrosi et al 2012) and adhesion (Baumgartner et al 2000).

This paper is organized as follows. Section 2 briefly introduces the TCAT framework. Section 3 describes the general mathematical formulation together with the constitutive equations. Section 4 briefly explains the computational model, and its numerical solution. Three examples of biological relevance are presented in section 5, and the conclusions follow in section 6.

The TCAT framework provides a rigorous yet flexible method for developing multiphase, continuum models at any scale of interest. An important feature of the procedure is that it explicitly defines larger scale variables in terms of smaller scale variables. When modeling transport in multiphase systems, the length scale of the model impacts the form and parameterization of the governing equations. At the microscale—the smallest scale at which the continuum hypothesis holds—a single (continuum) point contains a large number of molecules such that properties such as the density, temperature and pressure of a phase can all be defined. At the microscale, classical 'point' conservation equations and thermodynamic expressions are written. However, the domains of many problems of interest are too large, and the phase distributions are too complex to be modeled at the microscale only. The level of detail required to account for the geometric structure and the variability of variables at the microscale allows simulation of only very small domains. To overcome this challenge, many porous media models are formulated at a larger scale, called the macroscale—adequate for describing system behavior while filtering out the high-frequency spatial variability. The standard continuum mechanics approach to formulating these models is a direct approach wherein the conservation equations are written at the larger scale and a rational thermodynamic approach is employed to obtain closure relations. Although this approach can be mathematically consistent, the use of rational thermodynamics fails to retain a connection between larger scale variables and their microscale precursors (Maugin 1999, Jou et al 2001). Thus mathematical elegance is achieved typically at the price of inconsistent variable definitions and an inability to relate quantities at one scale to those at another scale. By averaging conservation and thermodynamic equations, TCAT avoids both of these pitfalls and leads to equations that are both thermodynamically and physically consistent.

The macroscale depends on the concept of the representative elementary volume (REV), an averaging volume that can be centered at each point in the system and which is large enough to include all phases present such that averages are independent of the REV size. The volume must also be sufficiently small so that quantities such as gradients are meaningful. TCAT consistently transforms microscale conservation and thermodynamic equations to the macroscale and converts averages of microscale derivatives into derivatives of macroscale average quantities. The description of a multiphase system must include dynamic conservation and thermodynamic equations for all phases, interfaces (where two phases meet), common curves (where three interfaces meet) and common points (where four common curves meet). Averaging theorems transform equations describing processes in these entities from the microscale to the macroscale (Gray et al 1993). Thus a physically complete system description is obtained that differs from the one obtained by only considering dynamic equations for phases with jump conditions at interfaces.

To close the conservation equations—which contain additional terms due to averaging—new model parameters and constitutive relations must be specified. While most methods use simplifications and the inclusion of approximate supplementary relations, TCAT employs averaged thermodynamic relations to guide closure of the system of equations. The microscale formalism chosen for averaging is the classical irreversible thermodynamics. This is adequate for the present model, but more complex formalisms can be employed (Gray and Miller 2005). The results are required to be consistent with the averaged entropy inequality that is also consistent with its microscale formulation.

The benefits of using a TCAT approach are as follows. Firstly, the model derivation proceeds systematically from known microscale relations to mathematically and physically consistent larger scale relations. This is accomplished by the use of averaging theorems. Secondly, the thermodynamic analysis is consistent between scales, in the definitions of variables at different scales and in satisfying the entropy inequality. The interscale consistency and explicit definition of variables are not achieved using a rational thermodynamic approach. Thirdly, relations may be obtained for the evolution of the spaces occupied by phases and of the interfacial area density. These relations are based on the averaging theorems. Although TCAT has heretofore been employed primarily in hydrology, it can impact tumor modeling in that the underlying physics and mathematics needed to describe tumors are related. Additionally, if hybrid tumor models are to be developed in the future, it is essential that the relation between the smaller scale variables and continuum variables be known. TCAT ensures that these relations are known.

The proposed computational model comprises the following phases: (i) the TCs, which partition into living cells (LTCs) and necrotic cells (NTCs); (ii) the healthy cells (HCs); (iii) the ECM; and (iv) the IF (figure 1).

Zoom In Zoom Out Reset image size

The ECM and IF pervade the whole computational domain, whereas the TCs and HCs are limited only to the subdomains with the tumor mass and healthy tissue, respectively. The ECM is modeled as a solid, while all other phases are fluids. The TCs become necrotic upon exposure to low nutrient concentrations or excessive mechanical pressure. The IF, transporting nutrients, is a mixture of water and biomolecules, as well as nutrients, oxygen and waste products. In the following mass and momentum conservation equations, α denotes an arbitrary phase, t the TCs, h the HCs, s the extra cellular matrix (ECM) and l the IF.

3.1. The governing equations

The governing equations are derived by averaging from the microscale to the macroscale and then using closure techniques to parameterize the resultant equations. The derivation is mathematically intensive such that providing it here in detail would distract us from the main thrust of this paper. The techniques have been employed for transport and for multiphase systems elsewhere (Gray and Miller 2009, Jackson et al 2009) and the procedure is the same for the current system, although the number of phases is different. An important feature of the approach is that the interphase contacts are explicitly accounted for.

The ECM is treated as a porous solid and porosity is denoted by ε, so that the volume fraction occupied by the ECM is ε^s = 1 − ε. The rest of the volume is occupied by the TCs (ε^t), the HCs (ε^h) and the IF (ε^l). Indeed, the sum of the volume fractions for all phases has to be unity:

The saturation degree of the phases is: S^α = ε^α/ε. Indeed, based on the definition of porosity ε and volume fraction ε^α in equation (1), it follows that

The mass balance equation for an arbitrary phase α based on the application of the averaging theorems is written as

where ε^α is the volume fraction, ρ^α is the density,  is the local velocity vector,  are the mass exchange terms accounting for transport of mass at the interface between the phases κ and α, and is the summation over all the phases exchanging mass at the interfaces with the phase α. However, if the interface is treated as massless, the transfer is to the adjacent phases, designated as κ. An arbitrary species i dispersed within the phase α has to satisfy mass conservation too, and therefore the following equation is derived by averaging:

where  identifies the mass fraction of the species i dispersed with the phase α, ε^α r^iα  is a reaction term that allows us to take into account the reactions between the species i and the other chemical species dispersed in the phase α and  is the diffusive velocity of the species i.

In particular, the mass conservation equation of the nutrient species i in the IF (phase l) reads

where it is assumed that no chemical reaction occurs within the phase and that the exchange of mass in the liquid is only with the tumor phase. Summing equation (5) over all species gives

where

Note that the mass exchange from the liquid to the tumor is actually to the living cell (LTC) portion of the tumor phase. The necrotic portion of the tumor is inert and does not exchange any nutrient with the IF. Also there is no need to make a distinction between the solvent part of the liquid phase and any of the dissolved species. All species are in the liquid phase. However, due to the relatively low concentrations of chemicals, the solvent phase is the dominant species and hence the global physical properties of the IF, such as density, intrinsic permeability and dynamic viscosity, are essentially those of the solvent.

The tumor phase t comprises a necrotic portion with mass fraction  and a growing phase with living cells whose mass fraction is . Thus the conservation equation for each fraction would be similar to equation (5). Assuming that there is no diffusion of either necrotic or living cells, and that there is no exchange of the necrotic cells with other phases, the mass conservation equation for the necrotic portion reads as

where ε^t r^Nt is the rate of death of TCs, or in other words the rate of generation of necrotic cells. In contrast to a mass exchange term between phases ( in equation (6) for instance), the reaction term ε^t r^Nt  is an intra-phase exchange term. The mass balance equation for the LTC is given as

where includes the exchange of nutrients and solvent from the IF to the tumor. Summation of these two equations yields an overall mass conservation equation for the tumor phase as

We can expand equation (8) by the use of the product rule and substitute it into equation (10) to obtain an alternative form of the necrotic species equation as

For the ECM and HC, the mass conservation equation becomes respectively

For the ECM and HC phases no mass exchange is expected with any other phase.

The momentum equation for the arbitrary phase α, including multiple species i, is

where is the body force, represents the momentum exchange from the κ to the α phase due to mass exchange of species i,  is the stress tensor and  is the interaction force between phase α and the adjacent interfaces. When the interface properties are negligible, this last term is simply the force interaction between adjacent phases. Given the characteristic time scales (hours and days) of the problem and the small difference in density between cells and aqueous solutions, inertial forces as well as the force due to mass exchange are neglected so that the momentum equation simplifies to

From TCAT, see appendix A, it can be shown that the stress tensor for a fluid phase is of the form , with p^α being the averaged fluid pressure and 1 the unit tensor, and that the momentum balance equation can be simplified to

where R^α is the resistance tensor.

3.2. The constitutive equations

No special assumption has been made yet for the constitutive behavior of the different phases, except for the fluid phases described by equation (16). In this paragraph, constitutive relations are explicitly presented for describing (i) the TC growth and (ii) the TC death, as a function of the nutrients' mass fraction and local mechanical stresses, for equations (6) and (8), respectively; (iii) the rate of nutrient consumption from the IF, in particular, to the LTCs, for equation (5); (iv) the diffusion of nutrients within the porous ECM, for equation (5); (v) the interaction force among the phases, for equation (15); (vi) the mechanical behavior of the ECM; and (vii) the differential pressure between the fluid phases.

The formulation presented in the above paragraph can be further simplified by assuming that the densities of the phases are constant and equal:

3.2.1. Tumor cell (TC) growth

This is regulated by a variety of nutrient species and intracellular signaling. However, without losing generality, in the present model a single nutrient is considered: oxygen. The case of multiple species can be easily obtained as a straightforward extension of the current formulation. TC growth is related to the exchange of nutrients between the IF and the living portion of the tumor. Therefore the mass exchange term in equation (6) represents tumor growth and, similarly to a part of the relevant equation in Preziosi and Vitale (2011a, b), takes the form

where the coefficient γ_growth^t  accounts for the nutrient uptake and the consumption of water needed for cell growth from the IF;  is the local mass fraction of the nutrient, a fundamental variable in the problem;  is a constant critical value below which cell growth is inhibited; and the constant  is the environmental mass fraction of the nutrient. Also, p^t denotes the TC pressure and its critical value p_crit^t  above which growth is inhibited. The Macaulay brackets  indicate the positive value of its argument. Note that, since the local nutrient mass fraction  within the tumor domain can be equal or smaller than , it derives that the non-negative part of the argument of the Macaulay brackets varies between 1 () and 0 (). Consequently, the growth rate for the viable TCs could at most be equal to γ_growth^t . Also in equation (18), H is the Heaviside function which is zero for p^t >p_crit^t  and is unity for p^t <p_crit^t . Note that  is the mass fraction of TCs that are necrotic and hence  is the volume fraction of viable TCs.

3.2.2. TC death

The rate of TC death in equation (8) can be described by the relation

where γ_necrosis^t  is the rate of cell death. All the other terms are similar to those presented in equation (18). However, the negative part of the argument of the Macaulay brackets  is considered. Also, p_necr^t  is the pressure above which the tumor stress has effect on the cell death rate, and δ_a is the additional necrosis induced by a pressure excess. Note that the mathematical form of equation (19) is very similar to equation (18) in that cell death is assumed to be solely regulated by insufficient concentration of nutrients (oxygen) and excessive mechanical pressure. No drugs or other pro-apoptotic molecules are used in the present model, but equation (19) can be readily modified to include also this contribution. Mathematically, a therapeutic agent or drug would be treated just as a 'nutrient'. Effects of cell membrane rupture and consequent transfer of liquid from the TC phase to the IF has not yet been included in the model. These aspects will certainly be included in future extensions of the current computational model. This will also be connected with the release of chemo-attractants and subsequent infiltration of macrophages. As such, these aspects can influence the local IFP.

3.2.3. The rate of nutrient consumption

As the tumor grows, nutrients are taken up from the IF so that the sink term in equation (5) takes the following form:

Nutrient consumption from IF is due to two contributions, namely (i) the growth of the TCs, as given by the first term within the square brackets in equation (20); (ii) the normal metabolism of the HCs, as presented in the second term. Indeed,  is related to the tumor growth, as discussed above; whereas the coefficient  relates to the normal cell metabolism. Being the nutrient mass fraction  in the tumor extra-cellular spaces always equal to or smaller than , the argument of the sine function varies between π /2  and 0. The part of consumption of oxygen related to the cells' metabolism depends on the oxygen availability and becomes zero when the mass fraction of oxygen is zero; this allows the values of the local mass fraction of oxygen to always be positive, since negative values have no physical meaning.

3.2.4. The diffusion of nutrients through the extracellular matrix

To approximate the diffusive flux in equation (5), Fick's law is used (). The effective diffusion coefficient of nutrients in the extracellular spaces is given as

where  is the diffusion coefficient in the unbound IF and δ is a constant coefficient greater than 1 which takes into account the tortuosity of the porous network. Actually the effective diffusion coefficient of oxygen does not have a linear dependence on the volume fraction of the IF, because it depends on the connectivity grade of the extracellular spaces. δ is a parameter that has to be calibrated experimentally.

3.2.5. The interaction force among the phases

R^α of equation (16) is the resistance tensor that accounts for the frictional interactions between phases. For example, porous medium flow of a single fluid encounters resistance to flow due to interaction of the fluid with the solid. If one has to model the flow at the microscale, a viscous stress tensor within the fluid phase would be employed. At the macroscale, the effects of the viscous interaction are accounted for as being related to the difference in velocities of the phases. The coefficient of proportionality is the resistance tensor. In multiphase flow, resistance tensors must be developed that account for the velocity differences between each pair of phases. Equation (14) contains the interaction vector  that arises between each pair of phases. In the full implementation of the TCAT analysis, the simplest result is that this vector is proportional to the velocity difference between the two indicated phases with the resistance tensor being the coefficient of proportionality. In the present version of the model, the interaction force  between the fluid phase α and the solid phase s (the ECM) is explicitly taken into account while the macroscopic effect of the interaction forces between the fluid phases , and is taken care of through the relative permeability k_rel^{α s} . The form of  is here assumed following the modeling of multiphase flow in porous media (Lewis and Schrefler 1998), that is to say

where k^{α s}  and μ^α are the intrinsic permeability tensor and the dynamic viscosity, respectively. Since there is no information available about this relative permeability which is a nonlinear function of the saturation and varies between 0 and 1, the following form is assumed

Equation (23) respects the constraint  and gives realistic results in agreement with the classical models present in the literature on porous media mechanics (Brooks and Corey 1964, Corey et al 1956, van Genuchten 1980). A more accurate determination of k_rel^{α s}  should derive from specific experiments or by the application of Lattice–Boltzmann modeling or analysis of micro-models. By introducing (22) in (16), the relative velocity of the phase α is derived as

The intrinsic permeability tensor k^ls  of the IF phase is constant and isotropic.

Experimental evidence confirms that cells would stay in contact with the ECM if the mechanical pressure gradients exerted over the cell phase are smaller than a critical value (Baumgartner et al 2000). For this reason, for the HCs and TCs the intrinsic permeability tensors (i.e. k^hs  and k^ts ) are isotropic but not constant, and are computed using the following equation:

This represents in mathematical terms the fact that if cells adhere firmly to the ECM, the phase permeability within the ECM is reduced. The minimum value of the permeability (set equal to ) eliminates the indeterminacy in the case , contained in the approach of Preziosi and Tosin (2009). This is an analogue in fluid dynamics to the stick-slip behavior in contact mechanics (Zavarise et al 1995).

3.2.6. The mechanical behavior of the ECM

The closure relation for the stress tensor acting on the ECM (sole solid phase) is

with  the effective stress tensor in the sense of porous media mechanics and the solid pressure p^s given as (Gray and Schrefler 2007)

where χ^α is the solid surface fraction in contact with the respective fluid phase, known as the Bishop parameter. This parameter is a function of the degree of saturation and is taken here to be equal to this last one (i.e. χ^α = S^α ). The Biot coefficient  is equal to 1 because of the incompressibility of the ECM. Indeed, this does not mean that the ECM cannot deform. The constitutive behavior of the solid phase is that of an elasto-visco-plastic solid in the large deformation regime (Zienkiewicz and Taylor 2000).

3.2.7. The differential pressure between the three fluid phases

The differential pressure between the fluid phases is a different concept from the interaction forces dealt with in section 3.2.5. In brief, the interaction forces are in play when there is flow. The different velocities of the different phases set up resistance forces between the phases. These are the interaction forces discussed above. Differential pressure, on the other hand, can exist even at equilibrium. It is not related to flow processes but is a statement that the pressures in adjacent phases can be different. In multiple fluid flow in porous media, this difference in pressures can be attributed to the curvature of the interface between fluid phases and to the surface tension. In the tumor system, the interfaces between phases are also capable of sustaining a jump in pressure between phases. In fact cells have surface tension which influences their growth and adhesion behavior (Dunlop et al 2011, Ambrosi et al 2012, Bidan et al 2012). At the microscale, the pressure difference between the cell phases and the fluid phases is equal to the interfacial tension, σ_c, multiplied by the interfacial curvature. After transformation to the macroscale, a macroscale measure is needed as a surrogate for the interfacial curvature. In porous media analyses, a surrogate for the pressure difference between fluid phases is proposed heuristically as a function of the fluid saturations (e.g., Brooks and Corey 1966, van Genuchten et al 1990). The cell pressure becomes very large when the available pore space is occupied by the cells, i.e. when S^l tends to zero. This behavior is depicted in figure 2. The following equation is proposed as a model for the pressure difference between the IF phase pressure p^l  and those of the cell phases p^t  and p^h

where σ_c and b are constants. The use of equation (28) to account at the macroscale for the curvature of the interface between the phases is an approximation that assumes that the distribution of the cells within the pore space does not impact on the pressure difference between the phases. This expression can be refined subsequently in light of experimental analysis.

Zoom In Zoom Out Reset image size

3.3. The final form of the governing equations

The primary variables of the model are: the tumor saturation S^t , the HC saturation S^h , the IF pressure p^l  and the nutrient mass fraction , together with the displacement of the solid phase (ECM) u^s . The time derivative of the latter is the ECM velocity v^s. After substituting for the explicit form of the constitutive equations, the final form of the governing equations is obtained. The mass balance equations of the ECM, TC, HC and IF are, respectively,

Summing equations (30)–(32), using the constraint equations on porosity and saturation, gives

The mass fraction of the necrotic cells is obtained from equation (11) as

The mass balance equation of the nutrient, using Fick's law to approximate the diffusive velocity () and assuming (17), is

Expanding equation (35) by the use of the product rule and substituting equation (6) gives an alternative form of the advection–diffusion equation of the nutrient species:

The linear momentum balance of the solid phase in a rate form (Schrefler 2002) is

where the interaction between the solid and fluids, inclusive of the cell populations, has been accounted for through the effective stress principle, i.e. equations (26) and (27).

Finally, for the solid phase the constitutive relationship between the effective stresses  and the elastic strains e_el^s , which is the difference between total strains e^s  and visco-plastic strains e_vp^s , reads as

where D_s  is the tangent matrix containing the mechanical properties of the solid skeleton. An extensive description of the elasto-visco-plastic model and of the mechanical aspects will be presented in a forthcoming paper.

The weak form of equations (30), (31), (33), (36) and (37) is obtained by means of the standard Galerkin procedure and is then discretized in space by means of the finite element method (Lewis and Schrefler 1998). Integration in the time domain is carried out with the generalized mid-point rule where an implicit procedure is used. Within each time step the equations are linearized by means of the Newton–Raphson method. For the FE discretization the primary variables are expressed in terms of their nodal values as

where , , ,  and  are vectors of nodal values of the primary variables at time instant t, and N_n, N_t, N_h, N_l and N_u are vectors of shape functions related to these variables.

For the solution of the resulting governing equations, a staggered scheme is adopted with iterations within each time step to preserve the coupled nature of the system. The convergence properties of such staggered schemes have been investigated by Turska et al (1994). In particular, for the iteration convergence within each time step a lower limit of Δt/h² has to be observed. Such a limit has also been found by Murthy et al (1989) for the Poisson equations and by Rank et al (1983) invoking the discrete maximum principle. The existence of this limit means that we cannot diminish at will the time step below a certain threshold without also decreasing the element size.

Three computational units are used in the staggered scheme: the first is for the nutrient mass fraction, the second is for computing S^t , S^h  and p^l  and the third is used to obtain the displacement vector u^s . Within each coupling iteration, equation (36) is solved for the mass fraction of the nutrient . Then the group of equations (30), (31) and (33) is solved in a fully coupled way for S^t, S^h and p^l. In this second computational unit, at each iteration i the approximate solution , , is used to update the mass fraction of the NTC , equation (34), the mass exchange term, equation (18), and the reaction term, equation (19). Once convergence is achieved for the second computational unit, the pressure in the cell phases (given by equation (28)) is used to compute the solid pressure, equation (27). The solid pressure is needed to solve the momentum balance equation (37). Once convergence is achieved within a time step the procedure can march forward.

Taking into account the chosen staggered scheme, the final system of equations can be expressed in a matrix form as follows, where some of the coupling terms have been placed in the source terms and are updated at each iteration to preserve the coupled nature of the problem.

with

where . The nonlinear coefficient matrices , and are given in appendix B.

The modular computational structure allows us to take into account more than one chemical species, simply adding a computational unit (equivalent to the first one used for the nutrient) for each of the additional chemical species considered.

The procedure has been implemented in the code CAST3M (http://www-cast3m.cea.fr) of the French Atomic Energy Commission taking advantage of previous work done on modeling concrete at early ages (Gawin et al 2006). There is a striking analogy between the two physical problems (concrete hydration and tumor growth) as far as the balance equations are concerned. In both we have one solid phase and immiscible fluid phases together with reactions and mass exchanges.

The computational framework above has been applied to solve three cases of practical interest: (i) the growth of a multicellular tumor spheroid (MTS) in vitro; (ii) the growth of an MTS in vivo; and (iii) the growth of a tumor along micro-vessels (tumor cord model). For all cases, the growth of the tumor mass, including the necrotic mass and living TCs; and the consumption of nutrient (oxygen) are analyzed over time. A direct comparison with the experimental data is presented for case (i). The ECM is assumed to be rigid for all three cases. This assumption will be relaxed in future studies. Results are presented in terms of volume fractions, ε^t, ε^h and ε^l, pressures, p^c  and p^l  and mass fraction of oxygen .

5.1. Growth of a multicellular tumor spheroid (MTS) in vitro

MTSs can be efficiently used to study the in vitro growth of tumors in the avascular stage. The tumor size can be easily measured experimentally using microscopy techniques and can be predicted quite accurately by analytical and computational methods. Here, the time evolution of an MTS is considered, assuming that the cellular mass is floating in a quiescent, cell culture medium. The geometry and boundary conditions of the problem are described in figure 3. Modeled as a half sphere imposing cylindrical symmetry, the MTS comprises three phases: (i) the LTC and NTC; (ii) the ECM; and (iii) the IF. At time t = 0 h, these phases coexist in the red area shown in figure 3, having a radius of 50 μm. Within this region, the initial volume fraction of the tumor cells (TCs) is set to 0.01, whereas the volume fraction of the ECM is set to 0.05 throughout the computational domain. Note that assuming a characteristic cell diameter of 10 μm, the initial number of tumor cells in the red area would be ∼ 10.

Zoom In Zoom Out Reset image size

The blue shell in figure 3—the cell culture medium surrounding the MTS—is the rest of the computational domain up to 1000 μm. These initial conditions are summarized in table 1. At the outer boundary (B1), the primary variables S^t ,  and p^l  are fixed with time (Dirichlet boundary conditions). At the symmetry boundaries B2, zero flux (Neumann boundary conditions) is imposed for all the phases. The atmospheric pressure is taken as the reference pressure. In this example, oxygen is the sole nutrient species, and its mass fraction is fixed to be  at B1 and throughout the computational domain at t = 0 h. The non-apoptotic cell death rate is calculated by equation (19), where the critical value of the oxygen mass fraction is given by , and the cell pressure above which the cell death rate increases is p_necr^t = 930 Pa . The necrotic regions are those where the mass fraction of necrotic cells  exceeds 0.5. All other governing parameters are listed in table 2.

Table 1. Initial conditions for an MTS in a medium.

	εs	εt	εh	pl
Red zone	0.05	0.01	0.00	0.00	7×10–6
Blue zone	0.05	0.00	0.00	0.00	7×10–6

Table 2. Input parameters used for simulating the first case.

Parameter	Symbol	Value	Unit
Density of the phases	ρ	1000	Kg m−3
Diffusion coefficient of oxygen in the IF		3.2×10–9	M2 s−1
Coefficient δ (equation (21))	δ	2.00	–
Intrinsic permeability for the IF phase	kls	1.8×10–15	M2
Intrinsic permeability for the tumor cell phase	kts	5×10–20	M2
Adhesion of tumor cells (to ECM)	at	1×106	N m−3
Growth coefficient of tumor cells (equation (18))	γgrowtht	0.016	–
Critical mass fraction of oxygen (equations (18), (20))		3×10–6	–
Critical pressure for cell growth (equations (18), (20))	ptcrit	1330	Pa
Necrosis coefficient (equation (19))	γnecrosist	0.016	–
The cell pressure above which necrosis increases (equation (19))	ptnecr	930	Pa
Pressure-dependent additional necrosis (equation (19))	δta	5×10–4	–
Consumption coefficient related to growth in equation (20)		4×10–4	–
Consumption coefficient related to metabolism in equation (20)		6×10–4	–
Coefficient σc in equation (28)	σc	532	Pa
Coefficient b in equation (28)	b	1	–

The volume fractions of the tumor cells ε^t  (TCs—solid line) and of the living tumor cells (LTCs—dashed line) with time are presented in figure 4(a). The radius for which the volume fraction ε^t  is zero gives the actual radius r^sph of the MTS. With time, the TC front moves outward and r^sph grows. The difference between the solid and dashed lines (TC–LTC) identifies the volume fraction of the NTCs. The LTC lines present a peak that moves outward with time, implying a continuous growth of the necrotic area within the MTS. Figure 4(a) clearly shows that r^sph grows from 50 μm (t = 0 h) to ∼ 400 μm at 360 h.

Zoom In Zoom Out Reset image size

The early evolution of the tumor mass is shown in more detail in figure 4(b). Starting from a 50 μm radius with S^t = 0.01, the tumor does not grow significantly in size within the first 50 h. The TC are rapidly dividing, increasing the volume fraction but not the size of the tumor mass. Thence, the tumor enlarges with a monotonic growth of r^sph. The tumor radius (r^sph(t)) is presented in figure 4(c) (solid line) along with the experimental data (open symbols). Notably, our prediction agrees well with three different MTS datasets (Yuhas et al 1977, Chignola et al 1995,  2000). The growth rate of the tumor is lower for the first 80 h. The numerical results are interpolated in figure 4(d) using the Gompertzian growth function

where r_∞ = 600 μm is the tumor radius r^sph at sufficiently large times (nominally t → ∞), and a and β are two constants derived from numerical data (a = 7.5 and β = 0.005 45 h⁻¹). The time t in equation (42) is measured in hours. For the necrotic core, a similar functional relationship is proposed here as

where is a constant, penalty term related to the thickness of the outer shell mostly comprised of viable cells (LTCs) which are still well nourished and oxygenated. The shell thickness depends on the cell line and nutrient availability (Mueller-Klieser et al 1986), but it is well accepted that at distances larger than 100–200 μm, nutrient diffusion is impaired. From our simulation, the shell thickness is 150 μm. Figure 4(d) also shows the necrotic core and the viable shell at three different times: necrotic cells are in the darker zone.

Note that the measured (apparent) volume V_app^TC  of the MTS could be very different from the effective volume V_eff^TC . Figure 5(a) shows these along with the effective volume of the LTC V_eff^TCL ; these are defined as

where , H is the Heaviside function which is zero for and unity for , and is the computational domain. The apparent volume also contains the IF, while the effective volume is comprised of TCs alone. Figure 5(a) shows that for short times, V_eff^TC and V_eff^TCL are equal as necrosis is initially negligible.

Zoom In Zoom Out Reset image size

The evolution of the oxygen mass fraction, the sole nutrient species considered here for cell proliferation and metabolism, is shown in figure 5(b). As the spheroid increases in size, gradients of oxygen concentration develop from the periphery, where the oxygen mass fraction is fixed at , to the center of the spheroid. Once the nutrient concentration in the center decreases below an imposed critical value (), cell necrosis commences. Note that at the boundary between the tumor and the surrounding cell culture medium, a significant change in the gradient of the mass fraction of oxygen is observed as a kink in the curves (see figure 5(b))—due to the lower effective diffusivity in the tumor. After a certain time, in the necrotic core the mass fraction of oxygen reaches a minimum value of about 1.5 × 10^–6, lower than the threshold critical value (). The oxygen concentration continues to fall below the threshold value until all cells are dead since in the necrotic region (here defined as the region where at least 50% of the cells are dead) the still living cells consume oxygen and slowly die. From figure 5(b), the necrotic core could also be identified as the portion of the MTS with a relatively homogeneous mass fraction of oxygen. The reasons: (i) necrotic cells do not consume oxygen, hence no nutrient gradient in the core; (ii) the local consumption of oxygen (regulated by equation (20)) decreases substantially and tends to zero with a reduction of the oxygen availability.

The pressure of the TC within the computational domain (computed using equation (28)) is plotted in figure 5(c). The maximum pressure of about 6.0 mm Hg (∼800 Pa) in the core of the MTS is lower than the critical pressure for cell death (p_necr^t = 930 Pa ). Thus oxygen deficiency is the sole cause of cell necrosis in the current example. Note that for relatively low saturations, the relationship between pressure and volume fraction is almost linear (see figure 2). Hence, the trends shown in figures 4(a) and 5(c) are similar. However, when the saturation level of the TC increases, the relationship with the pressure becomes nonlinear and so the peak pressure in figure 5(c) is more pronounced than the peak volume fraction ε^t (ε^t = S^tε) in figure 4(a). The IFP (p^l) is plotted in figure 5(d). Within the first 50 h, the TCs grow locally, while the overall external radius of the tumor mass stays constant at its original value (50 μm). As the IF is consumed by the TC, and the assumption (17) allows the volume balance to be satisfied locally the IF pressure gradient remains unaltered. Figure 5(d) shows that until 50 h the IF pressure gradient is zero so that no additional IF from the environment is needed (oxygen moves only by diffusion). After 50 h, the spheroid increases its radius; hence with tumor growth the IF must flow inward, per constraint equation (2). Therefore the IF pressure in the MTS core decreases. The intrinsic permeability of the IF phase is relatively high compared to that of the TC phase (see table 2). For this reason, the variations in pressure (figure 5(d)) are minimal but significant to explain that IF flows into the viable tumor shell during growth. Indeed, the IFP computed is slightly lower than in the surrounding tissue. This has to be ascribed to the lack of vasculature networks and lymphatic systems in the current model. The high IFP measured in tumors is mostly associated with the higher permeability of the fenestrated tumor endothelium and the lack of, or reduction in, lymphatic flow. Therefore the plasma permeating the tumor from the vascular compartment cannot be drained out efficiently through the dysfunctional lymphatic systems leading to the progressive liquid accumulation in the extracellular space and consequent pressure buildup (Jain and Stylianopoulos 2010). All this will be included in future extensions of the model also incorporating the vascular compartment and the lymphatic system. It should also be noted that, in the present computational model, the IFP depends strongly among other things on the pressure difference–saturation relationship of section 3.2.7 and possibly also on the deformation of the ECM. This aspect is currently under investigation.

5.2. MTS in vivo

In the second example, the tumor is growing within healthy tissue, which substitutes the cell culture medium in the previous case. Therefore, the initial configuration of the system comprises four phases: (i) the LTCs and NTCs; (ii) the host cells of the healthy tissue surrounding the tumor mass (HCs); (iii) the ECM; and (iv) the IF. The ECM and IF are distributed throughout the computational domain. The growing MTS pushes on the HCs as its radius increases.

Also, it is anticipated that the diffusion of nutrients towards the tumor mass would be reduced by the presence of the healthy tissue. As the tumor, the thin healthy tissue corona is assumed here not to be vascularized.

The geometry and boundary conditions of the problem are described in figure 6. The MTS is modeled considering a half sphere and imposing cylindrical symmetry. The red region contains the TCs with an initial radius of 30 μm (t = 0 h) and an initial volume fraction set to 0.45. The orange region is the healthy tissue extending till the outer boundary B1 of the computational domain of 150 μm. The volume fraction of the host cells in the healthy zone is initially homogeneous and set to 0.45 (t = 0 h). At B1, the primary variables S^t , S^h , and p^l  are prescribed and constant (Dirichlet boundary condition). At the boundaries B2, zero flux (Neumann boundary condition) is imposed for all phases and nutrients due to the radial symmetry. The atmospheric pressure is the reference pressure. As in the previous example, oxygen is the sole nutrient and its mass fraction is fixed at  on B1 and throughout the computational domain at t = 0 h. The chosen mass fraction of oxygen corresponds to the average of the dissolved oxygen in the plasma of a healthy individual. Although in this case the vasculature is not explicitly considered, the radius of the computational domain (here 150 μm) can be taken as an indicator of the vascularization grade of the host tissue: higher radii correspond to smaller vascularization and vice versa. Due to the lower reference environmental mass fraction of oxygen, the parameters γ_growth^t and , which govern growth and necrosis, respectively, and one coefficient of the oxygen sink term function (, see equation (20)) are different from the first example (see table 4). The initial conditions are listed in table 3, while the parameters of the healthy phase are given in table 4. All the other parameters are the same as in table 1.

Zoom In Zoom Out Reset image size

Table 3. Initial conditions for an MTS growing within healthy tissue.

	εs	εt	εh	pl
Red zone	0.05	0.45	0.00	0.00	4.2×10–6
Blue zone	0.05	0.00	0.45	0.00	4.2×10–6

Table 4. Additional input parameters for the second case.

Parameter	Symbol	Value	Unit
Intrinsic permeability for the host cell phase	khs	5×10–20	m2
Adhesion of the host cellsa (to ECM)	ah	1×106/1.5×106	N m−3
Growth coefficient of tumor cells (equation (18))	γgrowtht	0.0096	–
Necrosis coefficient (equation (19))	γnecrosist	0.0096	–
Consumption coefficient related to growth in equation (20)		2.4×10–4	–

^aIn the second case, the effect of cells adhesion is analyzed; then more than one value is used.

The adhesion of the cells to the ECM (a_t and a_h) has a more significant effect than in the first case, since there were no HCs surrounding the tumor mass. This example shows clearly that the relative cell adhesion plays a major role in affecting the overall tumor growth. The panels in figures 7(a) and (b) show the variation of the TC volume fractions over time for two different adhesion conditions, namely a_h = a_t (left column, figure 7(a)) and a_h = 1.5 a_t (right column, figure 7(b)). The times for the three panels are 1, 180 and 360 h. As for the first example, the variation of the volume fraction for the TCs is presented as a solid line, and the living portion of the tumor cells (LTCs) by a dashed line; the difference between the two gives the fraction of necrotic cells (NTCs). As expected, the overall tumor mass expands with time and a necrotic area appears in the core. Interestingly, the panels reveal a significant difference in the evolution of the volume fractions depending on the adhesion conditions. When a_h = a_t (figure 7(a)), the growing tumor mass completely displaces the HCs; whereas for a_h > a_t (figure 7(b)), the tumor spheroid during its growth pushes on the HCs and partially invades their domain. Interestingly, tumor invasion of the surrounding healthy tissue is controlled by the relative adhesion properties of the cell populations. Diffuse interface models and fourth order differential equations as used by Hawkins-Daarud et al (2012) are not needed for capturing the invasive behavior.

Zoom In Zoom Out Reset image size

Surprisingly, cell adhesion does not affect the overall volume size as clearly seen in figure 7(c) where the effective volumes of the tumor and the LTCs are plotted for the two adhesive conditions defined above. On the other hand, by comparing figures 7(c) and 5(a), the growth patterns of an MTS in a medium and an MTS in tissue appear quite dissimilar. This is mostly due to the presence of the adhesive host cell phase that contrasts the tumor growth and reduces the nutrient supply. In addition to the difference in growth pattern, a one to two orders of magnitude difference in effective tumor volume can also be observed. Hence the experimental results obtained in vitro are not indicative of the in vivo cases since the growth environments are very different.

The radius of the computational domain can be taken as an indicator of the vascularization grade of the host tissue, because at the boundary B1 the mass fraction is fixed at  (mass fraction of dissolved oxygen in the plasma of a healthy individual). We have solved the case with a_h = 1.5a_t for r_ext = 200, 250 and 300 μm (the values of the initial thickness of host cells, respectively, of 170, 220 and 270 μm) to evaluate the influence of the vascularization grade on the growth of tumor. The ratio between the effective tumor volume at 360 h and at t = 0 h has been plotted in figure 7(d) for different considered spherical domains: after 360 h the volume of the tumor is 18 times the initial volume for r_ext = 150 μm and 13 times the initial volume for r_ext = 200 μm; hence, if we increase r_ext the growth rate decreases.

5.3. Tumor growth along microvessels (tumor cord model)

In this last case, TCs grow in proximity of two otherwise healthy blood vessels that are the only source of oxygen. The presence of capillary vessels has an important impact on the tumor development and on its spatial configuration (Astanin and Preziosi 2009); this is confirmed in our application case where the progressive migration of TCs among adjacent vessels is also shown.

The system is comprised of four phases: (i) the LTCs and NTCs; (ii) the healthy tissue surrounding the tumor mass; (iii) the ECM; and (iv) the IF. The ECM and IF are distributed throughout the computational domain. The geometry and the boundary conditions of the problem are described in figure 8(b). We consider two straight blood vessels of 8 μm diameter. The TCs are initially located around one vessel only (see figure 8(a)). Two different separation distances between the vessels are considered: in the first simulation (S1) the distance is 80 μm; in the second (S2) the distance is 100 μm. Note that in these cases, a full 3D computational solution is required. The geometry has two planes of symmetry (i.e. the median horizontal plane and that passing through the two vessels, figure 8(b)); hence, only a quarter of the complete geometry is discretized. The FE mesh here is more complex than that of the previous cases. The parameters used are those of the second case, as treated in paragraph 5.2, with the exception of the volume fraction of the ECM, here set to 0.1. The mass fraction of the oxygen () is imposed as a boundary condition on the cylindrical surface of the two blood vessels. At the remaining bounding surfaces the flux of oxygen is zero. The fluxes of all the phases (l, h, t) are zero at the two symmetry planes and at the cylindrical surface of the two capillary vessels. For the remaining boundary the imposed conditions are shown in figure 8(b). The initial conditions are summarized in figure 8(a).

Zoom In Zoom Out Reset image size

The volume fractions at 7 and 15 days of the healthy cell phase HC and of the living tumor cells phase LTC are shown in figure 8(c) for case S1. The HCs are almost completely replaced by the tumor cells and after 15 days necrosis occurs in parts of the tumor which are more distant from the left blood vessel. Figure 8(c) shows also the oxygen mass fraction at 7 and 15 days for the same simulation. The strong decrease in the oxygen mass fraction, caused by the presence of the tumor, can be readily observed by comparing the areas populated by the abnormal and healthy cells.

In the second numerical simulation (S2), the distance between the two vessels is higher than in the case S1. The progressive expansion of the tumor mass from the left to the right vessel is therefore less. The results for the oxygen mass fraction are qualitatively similar to that of S1. Figure 9(a) shows the mass fraction of oxygen along a line passing through the two vessels, at different days. The effects of consumption of oxygen coming from the vessels are evident. Oxygen is here replenished only through the two blood vessels (two peaks). Figure 9(b) shows the tumor and the local mass fraction of oxygen (in the computational grid) at 15 days; clearly the higher values of the oxygen mass fraction are close to the two capillary vessels (). In S2, the tumor has not yet completely reached the right vessel after 15 days. In figure 9(c), the volume of the tumor after 20 days is represented for S2, and the necrotic area is clearly visible (only the finite elements in which the volume fraction of the tumor phase is higher than 0.01 are shown). In figure 9(d), the time evolution of the tumor volume is plotted for the two cases, S1 and S2. The plotted volume is that of the finite elements with a volume fraction of the tumor cells higher than 0.01. Note that initially there is no difference between the two cases because the growth is mainly influenced by the left vessel. After 10 days, the growth rate increases for the S1 case due to the additional nutrient supply coming from the right vessel.

Zoom In Zoom Out Reset image size

A tumor growth model has been developed based on multiphase porous media mechanics. The governing differential equations have been derived by means of the TCAT. These are mass balance equations for the different phases with the appropriate linear momentum balance equations. The equations have been discretized by means of the finite element method and a staggered procedure has been adopted for their solution. The lower limit of the ratio between the time step size and the square of the element size, necessary for a proper numerical behavior of staggered schemes and Poisson-type equations, has been determined by means of numerical tests.

The computational framework has been applied to three examples of practical interest, namely an MTS immersed in a cell culture medium, a tumor spheroid surrounded by healthy tissue and a tumor cord. Multiple phases have been considered in the model including (i) the LTC and NTC, (ii) ECM, (iii) IF and (iv) HC. For all cases, growth of the tumor mass, including the NTC and LTC areas, and the consumption of nutrient (oxygen) are analyzed over time within the whole computational domain.

For an MTS suspended in its culture medium, a direct comparison with three different experimental cases in the literature is presented. The agreement between the computational prediction of the tumor radius and the experimentally measured values is good. Also, the tumor growth follows the well-known Gompertzian growth pattern demonstrating again the accuracy of the computational model. Interestingly, the early development of the malignant mass is characterized by a rapid division of the tumor cells accompanied by an equally rapid increase in tumor cell volume saturation, whilst the overall tumor size stays almost constant. This was observed up to 50–60 h from the beginning. This early phase is then followed by a fast exponential growth (Gompertzian growth pattern). The model allows the volume of each individual phase to be calculated at each time.

In the second example, the MTS is surrounded by healthy tissue. The coexistence of two different cell populations (healthy and tumor) allows quantification of their relative adhesion to the ECM on tumor growth. In this respect, two different conditions are analyzed showing that when the HCs adhere less to the ECM, the tumor advancing front displaces uniformly the healthy tissue; in the opposite case the tumor cells infiltrate the healthy tissue at discrete points. Interestingly, this result has been achieved without involving diffuse interface models and fourth-order differential equations. The presence of the healthy tissue leads to an overall reduction in tumor growth, mostly due to the lower nutrient transport and geometrical confinement.

In the third example, the tumor cells along microvessels are predicted in a fully 3D geometry, with a clear delineation of the necrotic and living tumor regions. The progressive migration of tumor cells among adjacent vessels in search of additional sources of nutrients and oxygen is revealed. Also shown is that a larger distance between adjacent vessels needs a longer time for the tumor to grow, also demonstrating our model's capability to account for the vasculature.

The numerical accuracy and physical soundness of the computational model will increase the level of complexity that we can address in tumor biophysics—such as the contribution of the ECM stiffness, relative cell adhesion and IF pressure on the infiltration and development of malignant masses. Also, modeling the transport of therapeutic agents, in the form of individual drug molecules as well as nanoparticle, and angiogenic vascular growth, will be introduced in future extensions. A direct comparison of the predicted tumor behavior with experimental data derived from patients using clinically relevant imaging modalities should provide a validation of the presented approach. The modular structure of the framework allows straightforward inclusion of additional phases and nutrient types.

Faced with a continuously aging world population and the surge in cancer incidence, the approach presented here should engender novel therapeutic strategies and treatment optimization to improve the prognosis, outcome of intervention and quality of life for patients.

GS and BS acknowledge partial support from the Strategic Research Project 'Algorithms and Architectures for Computational Science and Engineering' (STPD08JA32—2008) of the University of Padova (Italy) and the partial support of Universita' Italo Francese within the Vinci Program. SS, WG and CM acknowledge partial support from the US National Science Foundation through grant number ATM-0941235 and the US Department of Energy through grant number DE-SC0002163. PD and MF acknowledge partial support from the NIH/NCI through grant numbers U54CA143837 and U54CA151668. MF acknowledges support through the Ernest Cockrell Jr Distinguished Endowed Chair.

The general conservation of momentum equation (14) will be denoted for the fluid phase using the letter f as a qualifier. In the paper, the f will be specified as either l, h or t.

where  is the body force,  represents the momentum exchange from the κ to the f phase due to mass exchange of species i,  is the interaction force between the phase f and the adjacent interfaces, and  is the stress tensor. If the inertial terms are considered to be negligible, as is the case for slow flow in a porous medium, the first two terms in equation (A.1) can be neglected. Additionally, the momentum exchange due to mass transfer  may also be considered small since this term is of the same order of magnitude as the inertial terms. Thus the momentum equation simplifies to

The TCAT method of closure involves arranging terms in the entropy inequality into force–flux pairs. At equilibrium each member of the force–flux pair will be zero. This equilibrium constraint guides closure of the conservation system for near equilibrium situations. In the case here where the flows are slow, the near-equilibrium state assumption is appropriate. Based on the TCAT procedure, the elements of the entropy inequality relating to flow velocity that arise in the entropy inequality are

In this equation,  is the macroscale temperature of the f phase,  is the gravitational potential,  is the chemical potential, p^f is the fluid pressure,  is the velocity of the solid phase and  is the rate of strain tensor of the phase f (). All of these quantities are macroscale averages.

Consider the variability in volume fraction of the f phase to be small. For this situation, . Additionally, consider an isothermal case such that the Gibbs–Duhem equation provides . Application of these two conditions to equation (A.3) reduces it to

This equation contains two independent force–flux products. The stipulation that both elements of each product pair must be zero at equilibrium and the requirement that the grouping of terms must be non-negative suggest the linear relations

and

In the first relation, R^f  is a symmetric, positive and semi-definite tensor accounting for the resistance to flow. In the second relation, A^f  is the fourth-order tensor that accounts for the dependence of the stress tensor on the rate of strain. At the macroscale for slow flow, this tensor is taken to be zero such that

is the resulting form of the stress tensor. We note that this does not imply that the fluid is inviscid. The effects of viscosity are accounted for at the macroscale by the momentum exchange term .

Substitution of the closure relations equation (A.5) and (A.7) into equation (A.2) provides the momentum equation in the form

Typically, this relation is expressed as

where  is called the hydraulic conductivity.

The hydraulic conductivity depends on the properties of both the flowing fluid and the solid porous material. For an isotropic medium, K^f = K^f 1 . The morphology and topology of the solid media are important in determining the hydraulic conductivity of the cellular solid phases. The conductivity is influenced by the cell size distribution, shape of the cells, tortuosity of passages, specific surface area and porosity (the sum of the fluid volume fractions). It also depends on the density and viscosity of the fluid. Neglecting gravity in equation (A.8) yields equation (16).

In the following equations K^s is the Bulk modulus of the solid skeleton and is the swelling strain rate depending on the solid pressure (equation (27)).